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A method for calculating the self-induced motion of a vortex sheet using discrete 
vortex elements is presented. Vortex panels and vortex filaments are used to simulate 
two-dimensional and axisymmetric vortex sheet roll-up. A straightforward application 
using vortex elements to simulate the motion of a disk of vorticity with an elliptic 
circulation distribution yields unsatisfactory results where the vortex elements move in 
a chaotic manner. The difficulty is assumed to be due to the inability of a finite number 
of discrete vortex elements to model the singularity at the sheet edge and due to large 
velocity calculation errors which result from uneven sheet stretching. 

A model of the inner portion of the spiral is introduced to eliminate the difficulty 
with the sheet edge singularity. The model replaces the outermost portion of the sheet 
with a single vortex of equivalent circulation and a number of higher order terms which 
account for the asymmetry of the spiral. The effect of the uneven sheet stretching is 
minimized by periodic rediscretization of the sheet. 

The resulting discrete vortex model is applied to both two-dimensional and axisym- 
metric sheets. The two-dimensional roll-up is compared to the solution for a semi-infinite 
sheet with good results. The axisymmetric cases show a smooth roll-up with none of 
the chaotic behavior seen in the straightforward approach. 

‘Adapted from the S.M. Thesis by John P. Kantelis, "Calculations of Axisymmetric Vortex Sheet 
Roll-Up Using a Panel and a Filament Model*, M.I.T., February, 1986. 
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Chapter 1 

Introduction 


A surface of discontinuity in the tangential component of the velocity field is called 
a vortex sheet. This three-dimensional surface, or sheet, can be thought of as a slip 
layer which is composed of vortex lines. A requirement for the existence of a vortex 
sheet is that the effect of viscous diffusion be very small, i.e., that the Reynolds number 
is very large. In the ideal case where the effect of viscosity is reduced to zero, the sheet 
is taken to be of zero thickness. In a real fluid the effect of viscosity is to immediately 
diffuse the vortex sheet into a shear layer of finite thickness, the rate of growth of the 
layer thickness being dependent on the fluid viscosity. 

In real fluids, vortex layers of finite thickness may form whenever a body with a 
sharp trailing edge such as a wing, or other lifting surface, moves through the fluid. 
This is the result of the generation of vorticity along the surface of the body ahead of 
the trailing edge due to viscosity. In the ideal fluid model, this condition is simulated 
by imposing the Kutta condition artificially at the sharp edge and the result is a vortex 
sheet of zero thickness. In the ideal case, it is also convenient to think of vortex sheets 
as being formed by moving a thin object through a fluid and then dissolving the object. 
At that instant there exists a vortex sheet in the shape of the object. An example of 
this is found in a paper by Taylor [22]. 

The most common examples of vortex sheets (or finite shear layers) are those which 
are shed into the wake at the trailing edge of an aircraft wing or helicopter rotor. The 
ensueing motion of these sheets is characterized by a rapid rolling up of the edges of the 
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sheet into concentrated vortex structures. Knowledge of the behavior of these sheets 
is desirable since they can have important influences on the behavior of the aircraft 
itself, particularly rotary-wing aircraft, as well as on nearby aircraft. The introduction 
of very heavy aircraft in recent years has renewed interest in the study of vortex sheet 
roll-up. These aircraft can produce very powerful rolled-up wake regions known as wing 
tip vortices which can be very dangerous to smaller aircraft should they happen to fly 
through them. Helicopter rotor performance is also greatly effected by these vortical 
regions since the loading on a helicopter rotor is strongly coupled with the motion of 
the rolled-up vortices in its own wake. 

This work is directed at developing a reliable and accurate method for calculating the 
motion of vortex sheets, particularly the motion and development of the sheet edge roll- 
up. Two-dimensional and axisymmetric sheets are studied with circulation distributions 
similar to that of the wake shed by an elliptically loaded wing or a helicopter rotor. The 
sheets are modelled using a discrete vortex element simulation where the continuous 
structure of the sheet is represented as a series of individual flat vortex panels or vortex 
filaments. The elements are followed in time in a Lagrangian manner with the velocities 
calculated numerically. 


1.1 Vortex Sheet Modelling 

Since the wake of an aircraft wing or helicopter rotor is characterized by very high 
local shear, it is a convenient and accurate approximation to model it as a vortex sheet 
immersed in a potential flow. Helmholtz’ theorems require that regions of vorticity in an 
inviscid fluid convect through the fluid with the local flow velocity. In an ideal fluid, any 
fluid element retains its vorticity as it moves. This is a very useful property since only 
the fluid particles having vorticity, namely the vortex sheet, need to be tracked as the 
flow evolves in time. This approach is commonly refered to as a Lagrangian description 
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of the flowfield since individual fluid particles of fixed identity are being tracked. A 
significant advantage computationally is realized here since typically in these types of 
flows only a very small region of the flow is of non-zero vorticity. 

The exact representation of the time evolution of a vortex sheet is given by a difficult 
singular integro-differential equation. For the general case the solution of this equation 
is not practiced and a numerical approach is taken. Typically, the continuous sheet is 
divided into a finite number of discrete elements, e.g., vortex filaments or flat panels, 
and each of these elements is tracked individually based on its local velocity obtained 
by the Biot-Savart law or some other method. 


1.2 Previous Work 

The first attempt at a flow computation for the roll-up of a vortex sheet using 
these techniques was performed by Rosenhead [18] in 1931. Rosenhead modelled a two- 
dimensional vortex sheet using a number of point vortices whose motion was tracked 
by computing the mutual induced velocities of all of the other vortices. It is perhaps 
fortuitous that Rosenhead had no computer available to do the arithmetic, forcing 
the work to be done by hand because it appears now that the fairly good results he 
achieved were in part due to the artificial smoothing of the data introduced by the 
manual computations. 

Later effort by Birkoff and Fisher [3] using a more refined computer analysis with 
more vortices than that used by Rosenhead resulted in much poorer results. The reason 
for this is thought to rest in the question of the existence and uniqueness of this initial 
value problem. It seems that an increasing number of point vortices to model the sheet 
does not yield a consistent solution of the original problem. It is not even clear that the 
inviscid problem is well-posed, as described by Saffman and Baker [20]. They indicate 
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that it is possible that smooth solutions may not exist for more than a finite amount of 
time. 

Another difficulty in the analysis of vortex sheets is their inherent instability. A 
two-dimensional vortex sheet is known to be unstable to small disturbances, the classic 
Kelvin-Helmholtz instability. This complicates the matter of computations for the roll- 
up of vortex sheets in that there are other reasons, related to the numerical procedure, 
why the sheet may become unstable. This instability, however, may be suppressed by 
artificially smoothing the sheet between time steps. Moore [12] shows this for the case 
of a circular vortex sheet. Smoothing has been applied indirectly in a number of vortex 
methods in the form of rediscretization of the sheet at regular intervals in time. This 
method was first used by Fink and Soh [7] , [8] where their point-vortex model of the 
two-dimensional vortex sheet is rediscretized at every time step to ensure that the point- 
vortices are equally spaced. The smoothing effect was a secondary consequence of the 
original motivation of Fink and Soh to rediscretize the sheet to reduce the truncation 
error in their velocity calculation. A numerical solution technique presented by Pullin 
[16] avoids the stability issue by computing the roll-up of a semi-infinite parabolically 
loaded two-dimensional sheet as a similarity solution. Pullin writes the governing equa- 
tion in a pseudo-time dependent similarity form which he solves iteratively for the sheet 
shape. The inner region of the resulting spiral agrees well with Kaden’s similarity solu- 
tion [10]. Pullin was also able to obtain the value of the arbitrary parameter inherent 
in Kaden’s solution. 

Additional two-dimensional models of vortex sheets have been presented including 
point-vortex models, vortex blob, cloud-in-cell, and panel methods. The vortex blob 
method [5] eliminates the singularity in the velocity field as the point-vortex is ap- 
proached by smearing out the vortex into a small finite region. This method is effective 
although somewhat arbitrary. The cloud-in-cell method [l] computes the velocities at 
each of the point vortices by solving the Poisson equation for the streamfunction and 
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then differentiating to obtain the velocity. The principle advantage of this method is its 
greatly increased computational efficiency, allowing for a greater number of vortices to 
simulate the flow. A difficulty with the method is the appearence of small scale struc- 
ture in the solution which appears to be grid dependent. The cloud-in-cell method has 
additional difficulties in axisymmetric flow where the self-induced velocities due to cur- 
vature are strongly grid dependent (Roberts and Murman [17]). Panel methods offer an 
advantage over point-vortex methods in that they yield a more accurate calculation of 
the induced velocity. Hoeijmakers and Vaatstra [9] presented a two-dimensional panel 
method which uses a splining technique for recomputing the sheet geometry at each 
time step. 


1.3 Motivation for this Study 

The purpose of this work is to extend the analysis of two-dimensional vortex sheets 
performed by authors such as Fink and Soh and Hoeijmakers into axisymmetric vortex 
sheets modelled with circular vortex panels and to compare the results with those ob- 
tained using circular vortex filaments. Axisymmetric sheets, being three-dimensional, 
have the added complexity of having additional self-induced velocity components. This 
arises from the well known fact that a curved vortex filament induces a velocity on itself. 
A combination of the methods developed for the two-dimensional sheets are employed 
and adapted for use on three-dimensional sheets. 

The validity of the method is demonstrated by computing the roll-up of a two- 
dimensional vortex sheet similar to the one shed by an elliptically loaded wing and 
comparing the results to the similarity solution for the semi-infinite sheet given by 
Kaden. 

The time evolution of two general configurations of axisymmetric sheet are studied, 
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the initially flat vortex disk resulting from a translating solid disk which suddenly dis- 
appears (i.e., Taylor’s problem), and a disk of vorticity more closely resembling that 
which would be produced by a helicopter rotor in hover. The translating disk is a good 
starting point for a study of this type since it is the axisymmetric analogue of the two- 
dimensional wake shed by an elliptically loaded wing. This later flow is certainly one 
of the most widely studied configurations and as such has a large assortment of results 
for comparison from previous studies in two-dimensions. The helicopter rotor case has 
more practical value where severed of the disks may be stacked on top of each other to 
simulate an evolving helicopter wake. 
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Chapter 2 

Problem Statement and Discrete Model 

2.1 Introduction 

The problem addressed here is that of an infinitely thin shear layer in an inviscid 
and otherwise irrotational fluid. The shear layer, hereafter refered to as a vortex sheet, 
convects through the fluid with its own induced velocity. Flows of this type are governed 
by Helmholtz’ vortex laws which state that the vorticity in the flow is conserved in time 
and that the vorticity is associated with material particles. That is, any fluid particle 
which had a certain amount of vorticity at the start of the flow will have the same 
vorticity for all time and any particle which was irrotational at the start will remain 
irrotational. Furthermore, it can be shown that the velocity field is uniquely determined 
by the distribution of vorticity. 

This is used to advantage here in that in order to follow the motion of the entire 
flow, only the motion of those particles which are rotational need be followed. This is 
very useful in that the rotational part of the flow is confined to a sheet. So by tracking 
only the convection of vorticity, the entire flowfield may be reconstructed. This type of 
analysis of motion is generally refered to as a Lagrangian analysis. 

The specific problem of interest is the motion of an axisymmetric vortex sheet, or 
disk. As shown by Batchelor [2], the velocity field, V(x ), is given by 

< 21 > 
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where s = x-x,f' is the strength of the sheet at the point x , and cr extends over the 
entire sheet. The resulting motion is then given by the solution of the integro-differential 
equation formed from Equation 2.1 with V(x) = dx/dt. 

This study will examine the motion of two axisymmetric sheets in particular, a 
sheet with a distribution of vorticity similar to that in the wake of an elliptically loaded 
wing and a distribution more typical of that in the wake of a helicopter rotor in hover. 
The initial motion of each of these cases is approximately known from the results of 
Kaden’s similarity solution [10] for a two-dimensional semi-infinite vortex sheet. This 
sheet has a vorticity distribution consistent with a parabolic loading which very strongly 
resembles the cases considered here near the edge. Therefore, in the initial time period 
the behavior should be similar. Kaden’s solution shows that the sheet immediately rolls 
up into an infinitely long logarithmic spiral at the tip singularity. It will be seen that 
this singularity is a cause of many problems in this type of analysis, however, it must 
be dealt with rationally since this feature is the driving force of the phenomena. 

The complexity of Equation 2.1 makes it unrealistic to expect an analytic solution 
for the vortex sheet roll-up, therefore a numerical approach is taken. The standard 
numerical method is to follow the motion of only a finite number of points along the 
sheet, each of which is assumed to represent the motion of an associated small piece 
of sheet, or, sheet segment. The problem is then to accurately calculate the induced 
velocities at these points, i.e., that which would result from evaluating Equation 2.1 
over each of the sheet segments. Due to the singular nature of Equation 2.1 and due 
to the curvature of the sheet the self-induced component of a segment on itself must 
be treated with some care. Given the induced velocities the points are tracked as they 
convect through the fluid. The two general approaches to sheet discretization have been 
to model the segments using vortex panels or vortex filaments. This will be done for 
both two and three-dimensional flows. 
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The vortex panel method involves discretizing the sheet into a finite number of seg- 
ments and then modelling each segment as a curved or straight panel of some predefined 
shape and strength distribution. The velocity field induced by an axisymmetric vortex 
panel is given by Sugioka [21] and is used here in a comparison with a filament model. 
The vortex filament method involves a further simplification in that all of the circulation 
contained in each segment is concentrated into a single vortex filament located at the 
center of the segment. The self-induced velocity of the segment is then taken to be that 
of the vortex filament. 

At luge distances from the sheet the induced velocities predicted by the panel and 
the filament models agree; however, as the sheet is approached, the flow field predicted 
by the panel model is smoother and more accurate. Both models display a large inac- 
curacy in the induced velocity field at close range due to the inherent ” graininess” of 
the discrete models. Whereas the vortex sheet is a continuous structure, the models 
considered here are discrete in that a singularity in the velocity field is present at each 
vortex or panel edge. 

Another method is to represent the circular vortex elements with sections of straight 
vortex filaments arranged into polygons. The velocity of each of the straight pieces is 
determined by evaluating the Biot-Savart integral over all of the remaining pieces. This 
is commonly known as the cut-off method. The advantage over axisymmetric panels or 
filaments is that the polygons will allow the sheet to develop in a fully three dimensional 
manner. The finite length of the elements introduces an effective core size with the 
representation of a vortex ring by a polygon. While not evaluated here, the relationship 
between a polygon and a ring is described in Appendix A. The polygon is the basis of 
free wake calculations that use vortex quadralaterals. The equivalence of these methods 
for axisymmetric flows will be considered in a later portion of the study. 
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2.2 Vortex Panel Model 


A detailed account of the induced velocity of a vortex panel is given by Sugioka [21] 
and is only outlined here. 

In the vortex panel model, the axisymmetric vortex sheet is subdivided into a number 
of flat axisymmetric panels or bands. Each of the panels has a total strength equal to the 
circulation of the sheet segment which it represents and is considered to be uniformly 
distributed over the panel. A description of a vortex panel and the coordinate system 
used in what follows is presented in Figure 2.1. The axisymmetric sheet geometry is 
described by a cylindrical coordinate system with r as the radial coordinate and z as the 
axial coordinate. All of the results presented here are for symmetry in the azimuthal 
direction. For the inner solution analysis, a panel centered coordinate system with an 
origin at the panel center is more convenient. The velocity component in the x or r 
direction is u and the component in the y or z direction is v. 

Sugioka obtained a uniformly valid expression for the induced velocity field of the 
panel by the method of matched asymptotic expansions. The solution approaches that 
of a circular vortex filament at large distances from the panel. Near the panel it closely 
resembles the two-dimensional solution for a vortex panel but includes (to first order) 
the effects of curvature. The velocity field is given as 

«(x, y) = «»n + tw “ «*o ( 2 - 2 ) 

v(x, y) = v in + v out - v io (2.3) 

where 

~ K (, * \ [ ... <?2 . n ( ti» - Qz W + Q 3 M 

“<” = 17 V 1 + 2R> r ^ Ql + V ~0T + “ ~qT) J 

k f x \ T Q 2 „ . . / ti7 — Q 3 tu + QsM 

Vi " = a (‘ + is) [- c ” 91o « qi - ( arctan -57- + arclan ~or)\ 
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and 


+2r[( w ~ Q») 1°ZQ* + ( w + Qs) log Qi - 4w 

t u — 0 < \ w + Q*, 

+2Q 4 arctan — — — - + 2Q4 arctan — — ] 
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Qi — Qs + 2aQs + to 2 
Qi = Qi ~ 2oQs + w 2 
Qs ~ x cos 6 + ysin0 
<?4 = ycos 0 — zsin0 

Qs = * 2 + y 2 


ri 

8/2 J 


r 1 = least distance from (z, y) to the panel center 
1*2 = greatest distance from (z, y) to the panel center 


where w is the half width of a panel, 2kw is the panel circulation or strength, R is the 
radius of the panel center, and K and E are the complete elliptic integrals of the first 
and second kinds respectively. 

The components «<,„< and v oxlt are the outer flow solution which is that of a single 
vortex ring. u,- n and v t - n are the inner solution components which are that for the two- 
dimensional panel including the effect of curvature and u t0 and v,- 0 are the components 
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which must be subtracted to complete the matching process. The primary advantage 
of the panel model over the filament model is that the resulting solution is uniformly 
valid in the entire domain including the panel itself (except near the edges) and is more 
accurate as the panel is approached. Additionally, the solution requires that w/R be 
much less than unity. This is violated as R -* 0 but since k -*• 0 as R — ► 0 this does 
not seem to pose a problem in the numerical method. 

The corresponding expressions for the induced velocity of a two-dimensional vortex 
panel are considerably simpler since the effect of curvature is not present. As given by 
Sugioka in the cartesian panel centered coordinate system of Figure 2.1, they are: 

u = — - sin 9 log — + 2 cos 0 (arctan — - -- + arctan 1 (2-4) 

4jt L Q l \ V4 V4 /J 

v — — — - cos 6 log ^ - 2 sin 6 (arctan — -pr — + arctan — ~t 1 (2.5) 

4ir L Q i \ _ Q 4 Sf4 / J 

where the Qi, k, and w are the same as for the axisymmetric case. These formulae are 
exact. 


2.3 Vortex Filament Model 

In the vortex filament model, the vorticity disk is modelled as a system of concentric 
circular vortex filaments, or vortex rings, which are initially evenly spaced along the disk. 
Each ring has a strength equal to the amount of circulation contained in the segment 
of the sheet which it represents. 

As is well known, a curved vortex filament possesses an infinite self-induced velocity 
at all points where its curvature is non-zero. In order to eliminate this infinite velocity 
it is necessary to consider the vortex rings as having a core of finite thickness with 
distributed vorticity. The self-induced velocity is now finite, but also dependent on 
the assumed distribution of vorticity within the core, which, unfortunately, is arbitrary. 
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However, since the vortex ring is intended to simulate the action of a vortex panel, a 
vorticity distribution may be assigned such that the two models become in some way 
equivalent. For example, as demonstrated by Sugioka [21], for equally spaced filaments 
a distribution of uniform vorticity within a circular core of diameter equal to 0.52 times 
the panel width (i.e., the distance between filaments) will yield equivalent kinetic energy 
and velocity for the two models. This is used as a criterion for determining the vortex 
ring structure. 

For a given core size, the associated velocity field of a vortex ring system may be 
determined by using the results presented in Lamb [11] for the induced velocity of a 
vortex ring. 

In general, the stream function at a point (r, z) due to a circular vortex filament is 
given by 

»('.*) = ~(r, + n)[«M - EMI ( 2 . 6 ) 

where: 


9 = Stokes stream function 
k = filament strength 
ri = least distance from (r, z ) to the ring 
rj = greatest distance from (r, z) to the ring 
<j> = (ri - r 2 )/(ri + r 2 ) 

K = complete elliptic integral of the first kind 
E = complete elliptic integral of the second kind 


The velocity components in the axial and radial directions are then 


v 


1 ££ 

r dr 


(2.7) 
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ia* 

u ~ r dz 


( 2 . 8 ) 


These are also applicable to the case of a vortex ring of finite core when the point (r, z) 
is not in or near the core. For (r, z) inside the core and at the center, the self-induced 
velocity of the single ring is 

8 R r 


v= — llog^-jj (2-9) 

u = 0 (2.10) 

where R is the radius of the ring, k is its strength, and c is the radius of the core. 

The two-dimensional form of the vortex filament is the familiar line or point vortex. 
In this case the velocity is given very simply as 

K 


U azimuthal — 


2 nr 


** radial — 0 


( 2 . 11 ) 

( 2 . 12 ) 


where k is the strength of the filament and r is the distance from the filament. In the 
cartesian coordinate system of Figure 2.1 the (x,y) velocity components are 


u = -k-t 


x 

V = K-r 


(2.13) 

(2.14) 


Then for either the panel or the filament model, the induced velocity at any position 
t is the summation of the contributions from each of the other elements and its self- 
induced velocity. The equations of motion for each ring or panel i become 

P-15) 

M J = 1 

(2- 16 ) 

M 3=1 

where for axisymmetric panels: 
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Uij — Equation 2.2 
Vi t j — Equation 2.3 

or for axisymmetric filaments: 


Ujj- = Equation 2.8 

when i ^ j 

= 0 

when i = j 

= Equation 2.7 

when * ^ j 

= Equation 2.9 

when t = j 


and for two-dimensional panels: 


Uij = Equation 2.4 


Vij = Equation 2.5 


or for two-dimensional filaments: 


Uij = Equation 2.13 

when * j 

= 0 

when i = j 

Vij = Equation 2.14 

when t ^ j 

= 0 

when i = j 


The motion of the sheet is given by the solution of the system of non-linear ordinary 
differential equations 2.15 and 2.16. 
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Figure 2.1: Vortex Panel and Coordinate Systems 
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Chapter 3 

Numerical Solution of the Evolution Equations 

The system of ordinary differential equations 2.15 and 2.16 governing the motion of 
the vortex filaments may be solved with a variety of numerical procedures. Considered 
here are the Euler method and a fourth-order Runge-Kutta method. Although generally 
more accurate, the Runge-Kutta method is much slower computationally then the Euler 
method, therefore, both methods will be evaluated to determine which has the best 
overall performance for this application. 


3.1 Euler Time Step Procedure 

The Euler method results from expressing the system 2.15 and 2.16 as the first order 
approximation to a Taylor Series. 

r!* +1 = r? + A trj* (3.1) 

z? +1 = + A tij* (3.2) 

where the superscripts refer to the time step level, t is the vortex element index, At is 
the time step, and f,- and are given by Equations 2.15 and 2.16. This method is easy 
to program and therefore suited to initial concept exploration runs, but if suffers from 
lack of accuracy unless the time step is very small. 
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3.2 Runge-Kutta Time Step Procedure 


The Runge-Kutta method results from modelling the solution over a time step as 
a polynomial rather than as a linear function as in the Euler model and is therefore 
more accurate for a given step size. However, the method requires the evaluation of the 
right-hand-side of the differential equations four times within each time step rather than 
once as in the Euler method. Therefore, for the Runge-Kutta method to be preferable 
over the Euler, the Runge-Kutta method needs to achieve accuracy greater that the 
Euler method for time steps that are four times larger. 


The Runge-Kutta method applied to Equations 2.15 and 2.16 yield the following set 

of algebraic equations to be evaluated at each time step: 

r „ +1 = rf + to + 2^ + 2 h+k, (33) 

<+■ = «r + +*i±k (3.4) 

where 

k 0 = Atf,(r" , 2 ”) 
ki = Atf,(r" + fco/2,z” + /o/2) 
k 2 = A tf,(r" + ki/2, z? + fi/2) 
k 3 = A tfj(r? + k 2 , z" + 1 2 ) 

and 

l 0 = Atij(r” , Zf) 
h = Atii(r? + k 0 /2, z? + Iq/ 2) 
l 2 = Atii(r” + k\/2, z f + fi/2) 

/ 3 = Atii(r i + k 2 , Zi + / 2 ) 

where f,- and i,- are given by Equations 2.15 and 2.16. 
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The value of the time step, At, is determined such that the maximum distance 
moved by any vortex element is limited to a specified value, A Xnw»*. Within each time 
step, after the velocities of all elements on the sheet have been calculated, the time step 
is computed as At = Ax max jV max where V max is the maximum velocity on the sheet. 

The relative merits of the two time step methods will be evaluated in Chapter 4 
where the vortex filament model is tested in a straightforward roll-up calculation. The 
time step procedures are evaluated based on their ability to conserve certain integral 
invarients of the flow; fluid impulse and kinetic energy. 
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Chapter 4 

Straightforward Application of Discrete Model 


4.1 Results 

A straightforward application of the vortex filament model to a vortex disk where the 
vortex filaments retained their identity throughout the calculation yields disappointing, 
although not unexpected, results. Similar to results for two-dimensional cases, the 
vortex rings begin to roll-up as expected, but then cease to define a sheet as smooth . 
lines connecting the vortex rings begin to cross over each other. A typical example of this 
for the vorticity sheet generated by the translating disk described earlier and modelled 
by 20 rings is presented in Figure 4.1. This case was run with the Runge-Kutta time 
stepping procedure with At = 0.005. The diameter of each filament core was chosen to 
be 0.52 times the width of the section of disk which it represents as discussed in Section 
2.3. Note that the sheet initially moves as expected but quickly the motion becomes 
chaotic with the sheet crossing itself. As discussed for similar cases in two-dimensions 
by, for example, Saffman and Baker [20], increasing the number of vortex rings results 
in no improvement. 


4.2 Effect of Time Stepping Procedure 


To eliminate the possibility of the chaotic motion being merely the result of inaccu- 
rate time stepping, the integral quantities mentioned in Section 3.2 may be evaluated 
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and checked for invariance. These invariants, impulse and kinetic energy, should remain 
constant in time. 

The total fluid impulse, given by P/p = ir £>r 2 summed over all of the vortex rings, 
is presented for the duration of the simulation in Figure 4.2 for the Euler method and in 
Figure 4.3 for the Runge-Kutta method. Each method was run with step sizes ranging 
from At = 0.1 to At = 0.0001. As expected, the Runge-Kutta procedure is far more 
accurate than the Euler procedure. A step size of At = 0.05 yields a nearly constant 
impulse for the Runge-Kutta method while the Euler method does not achieve constant 
impulse for even the smallest step size used. 


The total kinetic energy of the flow is given by Lamb [11] as 


n ot 


S 1/udzdr 


(4.1) 


where a / is the vorticity, 4' is the streamfunction, and p is the fluid density. Since this 
flow is irrotational everywhere except in the interior of the vortex ring, the energy can 
be expressed simply as a summation over each of the rings 

N 


T = —xp J ^ iOJidzdr 


(4.2) 


where N is the number of vortex rings and the integration is performed only over the 
filament core sections. For each core section t, the streamfunction is written as the sum 
of the contributions of all of the vortex rings including that of ring t itself. 


N 

*< = E 

i=i 


(4.3) 


where 'bij is the streamfunction at ring * due to ring j. With this the kinetic energy 
becomes 


N 

T=-«pZ 

•=1 


J J ^i'iUidzdr + J J *i,jVidzdr 


(4.4) 
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where is the streamfunction at ring t due to ring t. The first term in the brackets 
is given by Lamb as 


J J Vi.iWidzdr = - 


K}Ri 


2jt 


' 8 Ri 7 

log — “ 7 
Ci 4, 


(4.5) 


where «,• is the strength of ring i, Ri is the radius, and c t - is the radius of the core section. 
The second term is evaluated numerically by assuming that given by Equation 2.6, 
is constant over the core section of the ring. 


The resulting kinetic energy for the straightforward case is shown in Figure 4.4 for 
the Euler method and in Figure 4.5 for the Runge-Kutta runs. This shows similar 
results although smaller time steps are needed to achieve constant energy. The kinetic 
energy is not constant for the Runge-Kutta method until the step size is reduced to 
about At = 0.01, while again it doesn’t become constant at all for the Euler method. 

These results indicate that the Runge-Kutta time step procedure is not the cause of 
the observed chaotic motion. 


4.3 Suspected Cause of Chaotic Motion 

The reason for this chaotic behavior and lack of smooth roll-up is not known with 
certainty. The two usual arguments are (1) Kelvin-Helmholtz instability and (2) velocity 
calculation error. 

For the two-dimensional vortex sheet it is shown by Moore [13] that the sheet is 
unstable to small disturbances of any wavelength unless the sheet is rapidly stretching. 
That the sheet is most unstable to shortest wavelengths suggests that the initial value 
problem of two-dimensional sheet roll-up is not well posed. However, Safiman [19] ex- 
plains in a heuristic argument that the extension of the sheet due to stretching also 
extends the disturbances and, hence, their wavelengths. Since the rate of growth of the 
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disturbance decreases with increasing wavelength, the stretching can have a stabilizing 
effect on the motion. Moore and Griffith-Jones [15] and Moore [14] confirm this argu- 
ment, quantitatively, by determining that if the sheet stretches fast enough, such that 
the sheet strength drops off faster than time -1 / 2 , then the sheet is stable. Moore then 
applies the result to the tightly wound portion of a two-dimensional roll-up spiral with 
a parabolic distribution of circulation and finds that the condition is met for stability. 

For parts of the sheet where rapid stretching does not occur Moore [12] shows that 
a discrete form of Kelvin-Helmholtz instability can ruin the calculations by amplifying 
the round-off and truncation errors inherent in the numerical procedures. 

The calculations presented for the straightforward case show that the calculations 
break down in the region of the spiral roll-up where the stretching is greatest and, hence, 
the sheet is most stable. In addition, the region where the sheet is still nearly flat shows 
no sign of Kelvin-Helmholtz instability. Therefore, it is assumed that the cause of the 
chaotic motion is not due to Kelvin-Helmholtz instability. 

Fink and Soh [7], [8] have demonstrated for the two-dimensional sheet that large 
errors are introduced into the velocity calculations as the vortices move away from the 
centers of their segments. This situation occurs very quickly at the outer edge of the 
sheet where the roll-up takes place. The high concentration of vorticity in this region 
induces very high rates of stretching and moves the vortices away from the centers 
of their segments. Fink and Soh suggest that this is the cause of the chaotic motion 
observed so quickly during the roll-up process. It is assumed that the same type of 
velocity errors are present in the three-dimensional sheet calculations and that this is 
the immediate cause of the error in the roll-up calculations presented in Figure 4.1. 

A second major source of error in the straightforward approach is in the inability 
of a finite number of discrete elements to model the more tightly wound regions of the 
spiral. As can be seen from Kaden’s similarity solution for the roll-up of a semi-infinite 
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two-dimensional vortex sheet with a parabolic circulation distribution, the inner part 
of the core assumes a logarithmic spiral shape: 

r = aO v (4.6) 

where r and 6 are polar coordinates with origin at the center of the spiral, a is a scaling 
factor, and v —* — 2/3 as 0 — ► oo. At any time after the start of the motion the sheet is 
infinitely long and is rolled into loops of rapidly decreasing spacing. It seems doubtful 
that any discrete model could successfully represent this region of the flow. 

In summary, the straightforward approach has two major areas of difficulty: 

1. Arbitrarily large velocity errors result from stretching. 

2. Tip region not modelled correctly by a finite number of discrete elements. 

These problems are discussed in the next sections. 
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Chapter 5 

Special Models 

5.1 Model of Singularity at Sheet Edge 

5.1.1 Motivation for Special Model 

A method for reducing the severity of the errors realized in modelling the tightly 
wound portion of the sheet edge in a discrete manner that has been used with some 
success by, for example, Hoeijmakers and Vaatstra [9] is to replace the inner portion of 
the spiral with a single filament of equal strength. As the inner layers of the sheet become 
more closely spaced, they are amalgamated into a single tip vortex ring such that the 
total strength and centroid of the system is preserved. This process is generally refered 
to as core dumping. Criteria for dumping that have been used include: restricting the 
spiral to a maximum number of loops, requiring a minimum number of vortices per 
loop, and amalgamating when the spiral layers move too closely to each other. 

Another motivation for replacing the outer portion of the sheet with a single ring is 
more fundamental. It is doubtful that a finite number of discrete elements would ever 
be able to accurately model the singularity in the flowfield at the tip. This singularity is 
responsible for the initiation of the roll-up at the disk edge. In fact, the roll-up may be 
considered to be the initial condition for the motion of the sheet. The infinite velocity 
singularity at the sheet edge is resolved by the sheet assuming the shape of a spiral of 
infinite arc length placed at the sheet edge at time t = 0 + . The discrete model discussed 


41 


here, i.e., a series of equally spaced panels or filaments, does not in any way treat this 
aspect of the flowfield. In fact it may be questioned how any discrete model of a finite 
vortex sheet could show the proper type of roll-up behavior. 

5.1.2 Inner Spiral Model 

The asymmetry of the spiral cross section is neglected when using only a single vortex 
filament to represent the inner part of the spiral. A more accurate representation of the 
sheet edge and inner spired may be obtained by replacing it with a numerical model of 
the Kaden [10] similarity solution for the roll-up of a parabolically loaded semi-infinite 
two-dimensional sheet. The asymmetric component included by the Kaden model has 
an effect on the flow at all points outside of the inner spiral and on the central vortex 
filament itself. Therefore, there are two additional velocity components to be modelled 
when using a single vortex filament to represent the circulation of the inner spiral, (1) 
that of the inner spired asymmetry on the outer flow, and (2) that of the inner spiral 
asymmetry on the central vortex filament. 

The inner spiral model is illustrated in Figure 5.1. The inner part of the spiral, i.e., 
the portion of the sheet contained within the circle of radius R a , is considered to be 
replaced with an equivalent flow system consisting of a single vortex filament of equal 
strength, k, and a number of higher order flow elements. The filament, located at the 
center of the spiral, accounts for all of the circulation in the circle R' a and is also used 
to track, the motion of the spiral center itself. The higher order elements model the 
asymmetry of the spiral and have no circulation. These act on both the sheet outside 
of the circle R' a and on the central vortex filament. 

The additional asymmetric velocity components are time dependent. As the roll- 
up progresses, the inner layers of the spiral move more closely together reducing the 
asymmetry and its effect on the flow. This may be seen by considering the sheet at 
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two different times, t\ and <2 where ti > ti as shown in Figure 5 . 2 . The solid lines 
indicate the part of the sheet that is modelled using vortex panels or filaments while 
the broken lines indicate the part that is modelled by the inner spiral model. Since 
the points on the sheet are material particles, the edge of the sheet is always identified 
by the same fluid particle. In Figure 5.2 this particle is indicated as point A. The 
nature of the roll-up requires that point A always remain at the same distance from the 
spiral center, although it revolves around it. This constant distance, R' a , is used as a 
reference parameter. Observe that the two inner spirals have the same diameter but that 
the configuration at time t = shows the spiral layers much closer together. This trend 
toward symmetry demonstrates the time dependency of the inner spiral model. As time 
increases and more turns are formed in the spiral, the contribution of the asymmetric 
components approaches zero. 

The flow components due to the asymmetry of the inner spiral may be found by using 
the results of Puffin [ 16 ] for the roll-up of a parabolically loaded two-dimensional semi- 
infinite sheet. Puffin studied the roll-up by transforming the unsteady-time dependent 
problem into an integro-differential form using a similarity parameter involving time 
and circulation. This eliminated the need to study the problem in an unsteady manner 
along with the associated instabilities. 

The initial condition for the semi-infinite problem is a flat sheet with a circulation 
distribution given as T = 2 <iy/x where x is the distance from the edge to a point on 
the sheet and a is an arbitrary scaling factor. Puffin describes the position of the sheet 
with the complex function zo(r,t) where t- is time and T is the circulation which also 
identifies a particular fluid particle. Based on a dimensioned analysis, the solution is 
written in the form 

*o(r,<) = (at) 2 / 3 n(A) (5.1) 

where fl(A) = £(A) + 177(A) is the non-dimensional sheet shape and A is the similarity 


43 


parameter, 

A= _J[_ 

(aH ) 1 / 3 

Pullin solves for the sheet shape function ft (A) numerically using an iterative method 
which yields a similarity solution for the roll-up. 

He then goes on to write Kaden’s similarity solution in a form consistent with his 
notation as 

n-n. = ar-V , «p[<( 5 ^ ? + «)] (5.3) 

where 

r = A" 3 (5.4) 

ft„ is the location of the spiral center, and a and e are arbitrary constants. The values 
of a and e are determined by matching the numerical results to Equation 5.3 at a point 
in the inner spiral. The values found by Pullin sure a = 0.124 and e = 2.69. 

Pullin’s results are used to find the asymmetric velocity components by numerically 
integrating over the sheet inside the circle R.' A . For example, the velocity at the central 
vortex filament is found directly by integrating over the inner spiral for both the x and y 
cartesian velocity components. The flow external to the circle is determined somewhat 
less directly. An integration is performed over the sheet to yield the streamfunction 
on the circle. From this is subtracted the streamfunction of a point vortex of strength 
equal to that of the sheet yielding the streamfunction of only the asymmetric part. The 
streamfunction in the flowfield is then given as the solution to V 2, t = 0 with the value 
on the circle as the boundary condition. 

Once these solutions have been found, it is necessary to express them in a rotating 
frame of reference since the integrations are performed with the point A on the negative 
y axis as shown in Figure 5.3, but its true orientation will vary with the rotation of 
the inner spiral during the roll-up. This is one aspect of the time dependency of the 
problem mentioned earlier. The other is that the sheet inside of point A becomes more 
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symmetric as the roll-up progresses. Therefore, the numerical integrations must be 
performed over a range of stages in the roll-up to determine the time dependency. 


The numerical integration on the inner spiral is performed as follows. The sheet 
inside the circle of radius R A is divided into a number of straight segments, as shown 
in Figure 5.3, which are modelled as individual point vortices for the purposes of this 
integration. Associated with each point vortex i is a streamfunction and velocity com- 
ponents given by 

. IT *» 

(5.5) 


T » K i , r 

*‘ = ~2 7 log iS 


vl =*- 

2xr 


(5.6) 

K = 0 (5.7) 

where the primes indicate dimensional quantities, k' { is the strength of the filament, R a 
is the radius of the circle enclosing the spiral, and r is the distance from the filament 
to the point where the integration is performed. The strength, k[, is obtained from 
Pullin’s similarity parameter as 

= r|, - 1-;, (5.8) 

where and rj are the values of the circulation at the end points of sheet segment «. 
Then from Equations 5.2 and 5.4 


«: = - ^ ,/j ) 


(5.9) 


which when substituted into Equations 5.5, 5.6, and 5.7 yields the contribution of seg- 
ment » to the non-dimensional streamfunction, 'F, 


r. -1 / 3 — r- -1 / 3 

*. = _2 = M W- 

* ( a 4 *) 1 / 3 2 jt g 

and to the non-dimensional velocity components 

d' r-i/s _ _-l/3 

v B = vl Ka - ^ 

9i 9 < 


(5.10) 


(a 4 !) 1 / 3 


2rrr 


(5.11) 


45 


V ri = 0 (5.12) 

The velocity components are rotated to yield the non-dimensional components in the 
(x,y) system 

r . -1 / 3 — r.'i / 3 «_ - 

(5.13) 


- 1/3 - 1/3 

T ii ~ r u yc-yi 

U i = i — 4 

27rr r 

- 1/3 - 1/3 _ 

_ 7 »a g « 

’ 2?rr r 

where r = r'/R' A and is expressed in Pullin’s non-dimensional coordinates as 


(5.14) 




(5.15) 


In this expression, the subscript » refers to the segment of the sheet and the subscript 
c refers the the location of the point at which the streamfunction or velocity is calculated. 
( i and rji are obtained from Equation 5.3 as 

& - c* = « r r 2/s c ° s (jj~2 + c ) ( 5 ie ) 

»?,-»7« = 0(T-r 2/3 sin(^2+«) (517) 

where and rj v are the positions of the spiral center and r,- is the value of r at the 
center of segment i. 


In the numerical integration, each sheet segment, t, is chosen such that it subtends 
an angle A $ along the sheet. The values of r,-, r,-,, and r,- a are then obtained by noting 
from Equation 5.3 that r is proportional to the angle 9, in fact, Ar = 2ira 2 A0. 

The radius of the circle, R ' a , is constant and is given by R.' a = i?(at) 2 / 3 where 
R is the non-dimensional radius at time t. In an actual roll-up calculation, this may 
be evaluated once at time to when the simulation begins where Rq is given by i?o = 
\/(£o — £«) 2 + (*?o — »7«) 2 - Co and rjo are the positions of the edge of the inner spiral 
sheet at time to which are obtained from Equations 5.16 and 5.17 with ro = 0.24, the 
value at the bottom of the spiral. 
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Equation 5.10 is used to calculate the contribution to the streamfunction at a point 
on the circle R' a due to one segment of the spiral and Equations 5.13 and 5.14 to 
calculate the velocity at the spired center. The total effect due to the spiral is obtained 
by integrating over the entire sheet inside the circle. However, since the sheet is infinitely 
long, the calculations are ended as soon as the solution converges. The entire numerical 
integration over the inner spiral is written 



r~r A 



2 j r 


logr 


oo -1/3 —1/3 

~ _ y' h T ii Vc ~ Vi 

2xr r 


r=r A 


oo ,-1/3 _ _ 1/3 

v= ^2 — 6 “ ’ 

t=t A 


2nr 


(5.18) 


(5.19) 

(5.20) 


The dependency of the problem on time becomes apparent at this point. The limits 
on the integration of Equations 5.18 through 5.20 are from the outer edge of the sheet, 
r = T A> to the center of the spiral, r — ► oo. As time advances and the roll-up progresses 
the value of rx changes since rx = a*t/T A and Tx is constant. To determine the effect 
of increasing time (and symmetry) the calculations of Equations 5.18 through 5.20 are 
repeated over a range of values of rx, from the smallest expected to a very large value 
to determine the asymptotic behavior. 


The resulting calculations for the non-dimensional velocity components at the central 
tip filament are found to very closely fit the forms 


u = -0.000332r^ 2 ' 383 

(5.21) 

v = -0.004885r^ 1543 

(5.22) 


where (z, y) is the rotating coordinate system fixed to the sheet center. The dimensional 
components in the usual cartesian coordinate system are obtained by rotating Equations 
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(5.23) 


5.21 and 5.22 through the angle a rot and multiplying by (a 4 t) l / 3 /R' A 


u = (u cos a ro t - v sin a rot ) 


(o 4 t) 1/3 

R\ 


I 

V 


(u sin a ro t + v cos oc r ot) 


(o 4 0 1/3 

R'a 


(5.24) 


The angle a ro t is dependent on time and is the angle through which the (x, y) coordinate 
system must be rotated in order to match the point A on the inner spiral to the point 
A on the actual sheet (i.e., the edge of the actual sheet). This is illustrated in Figure 
5.4. The sheet edge position A is determined during the roll-up calculations. 


The flowfield external to the circle is handled somewhat differently. Equation 5.18 
is used to compute the asymmetric component of the streamfunction on the boundary 
of circle R,' a , denoted as 9 boundary The flow exterior to the circle is then obtained as 
the solution to V 1 ® = 0 with ^boundary as the boundary condition. The solution may 
be written as 


, ^ , x sin nO - , . cos nO . „. 

* = £ B n (r A )-— + D n (r A )-— (5.25) 

r»=l L r 

where B n and D n are written to emphasize their dependence on time. As shown in 
Figure 5.4 6 is measured from the negative y axis in the rotating (x,y) coordinate 
system. S n (rx) and D n (rX) are given by 


r n f 2 * . - ~ 

B n — ~j( Jq ^ boundary Sill T 10 d$ 

r n ft* > - 

Dn = “ ^ boundary COS ndd9 

where V boundary is also dependent on r*. B n {jA) and Z? n (rx) are found to be fit closely 
by the forms 

B n = b n T m " 


Dn = d n T l 


n 


where 
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n 

b n 

m n 

d n 

In 

1 

0.005301 

-1.386 

0.0010070 

-2.373 

2 

0.001310 

-1.363 

0.0001550 

-2.282 

3 

0.000575 

-1.345 

0.0000657 

-2.409 

4 

0.000319 

-1.330 

0.0000312 

-2.361 

5 

0.000202 

-1.316 

0.0000157 

-2.173 

6 

0.000139 

-1.301 

0.0000102 

-2.221 

7 

0.000100 

-1.283 

0.0000115 

-2.807 


are the first seven sets of coefficients. 

Although the numerical computations yield decimal numbers for the time behavior 
of the spiral solution, in reality these should be rational numbers. Based on the results 
obtained, possible rational exponents are —4/3 for the m n and —7/3 for the 


The non-dimensional velocity components in the flow field are then given by 


v . = -^- + JL. 

9 dr 2 st 

V . = l J± 

T r d8 


(5.26) 


(5.27) 


where (f , $) is the coordinate system fixed to the center of the spiral as shown in Figure 
5.4. The dimensional components in the usual cartesian coordinate system are given by 


' d . . . .( a 4 *) 1 / 3 

u = (Vf cos 8 t - V- t sin 9 t ) f — 

k a 

, f a 4 f U/3 

v = (Vf sin 6, + Vj cos 0 t ) - — -7 — 

Ra 


(5.28) 


(5.29) 


where 0 t = 6 — x/2 + a ro t and is taken positive in the usual sense with 6 t = 0 on the 
positive x axis. 
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5.1.3 Initial Representation of the Sheet 

The Pullin solution may also be used as an initial condition so that the numerical 
simulation may be started at some time, t = to, where to > 0. The outer most portion 
of the sheet will have already rolled up into an infinitely long logarithmic spiral which 
can be approximated using the Pullin model described in the previous section. This 
approach eliminates the problem of attempting to model the initial motion of the tip 
singularity in a discrete manner with a series of vortex elements. 

The state of the sheet at time t = to may be obtained from Pullin ’s results as follows. 
The outer portion of the sheet will have assumed the configuration shown in Figure 5.5. 
At time to, point Ao has convected to position A on the curved sheet and the location of 
the spired center is at ao, bo relative to the original location of the sheet edge. Pullin’s 
results yield the location as ao = 0.308(ato) 2 ^ 3 and &o = 0.489(a<o) 2 ^ 3 where to is the 
initial time and a is a measure of the total sheet strength as given in Section 5.1.2 for 
the parabolically loaded sheet as T = 2 ay/x. 

The approach taken here is to start the simulation at time t = to with the outer 
portion of the sheet from the edge to point Ao, and its associated circulation, represented 
as a single vortex filament and higher order components centered at location ao, &o- The 
radius of the core section of this circular filament is determined by preserving the total 
kinetic energy of the system. Also, for simplicity, the inner portion of the sheet at 
time t = to, i.e., from point A to the central axis, is modelled as a flat sheet that 
is stretched uniformly from point A to point Ao. This stretching procedure does not 
precisely conserve the fluid impulse but for the small values of to used in the simulations 
the impulse is nevertheless correct typically to within less than 1%. 
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5.1.4 Tip Ring Amalgamation 


As the sheet rolls up and the individual layers of the spiral come close together 
there will be large errors in the induced velocity calculated at the vortex elements. 
As a concentrated vortex is already present at the center of the spiral from the initial 
condition it is natural to use core amalgamation to lump together the inner portion of 
the spiral as the layers move too close together. 

Core amalgamation may be accomplished in between time steps by requiring the 
before amalgamation and after amalgamation systems to be equivalent. The strength 
and location of the resulting single ring may be determined by conserving the circulation 
and centroid of vorticity of the group. 

It is also necessary for the self-induced velocity of the centroid of the group and that 
of the amalgamated single ring to be equal. Bernardinis, Graham, and Parker [6] have 
expressed doubts about the amalgamation process when applied to three-dimensional 
flows unless this velocity can be maintained through the amalgamation process. In 
their work with vortex rings they do not explicitly require the velocities to be equal 
but instead determine the core radius (and hence the self-induced velocity) based on 
the conservation of the volume of the rotational fluid of the group of vortices. The 
approach taken here for the determination of the core radius is to conserve the kinetic 
energy of the flow across the amalgamation by adjusting the core radius of the resulting 
tip ring. In this way the three primary invarients of the flow, circulation, impulse, and 
kinetic energy Me conserved throughout the simulation. The amalgamation process is 
summarized in Figure 5.6. 

Limiting the spiral portion of the sheet to a fixed number of loops is the amalga- 
mation criteria used here. Another method would be to check for a minimum distance 
between loops since this is the reason that amalgamation is used, but it is assumed that 
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simply counting loops provides a good enough approximation. 

Another reason for limiting the number of loops is in consideration of the total 
length of the sheet. As more loops are allowed to form around the tip vortex ring 
core the sheet becomes excessively long requiring very many discrete vortex elements 
for sufficient resolution. The number of vortex elements, in turn, is limited by the 
computational resources available. 


5.2 Rediscretization to Alleviate Velocity Errors 

5.2.1 Motivation for Rediscretization 

Fink and Soh [7], [8] show that the error in calculating the induced velocity at a 
point on a two-dimensional vortex sheet can become arbitrarily large if the point is not 
in the center of the segment. This is the result of the neglect of a logarithmic term 
in the solution of the Cauchy principle value integral representing the velocity. When 
the vortex is centered in the segment, the logarithmic term vanishes. As the vortex 
approaches the edges of the segment, however, the contribution of the logarithmic term 
becomes large and cannot be neglected. Even if the sheet is initially discretized into 
equal length segments, the ensueing stretching and deformation of the sheet results in 
their movement off center. 

Fink -and Soh eliminate this error by rediscretizing the vortex sheet at every time 
step and achieve good results. This involves constructing a new smooth sheet geometry 
through the convected positions of the point vortices using a curve fitting procedure 
and then counting off equally spaced intervals along the new sheet to form the (equal 
sized) sheet segments for the next time step. The point vortices are then put at the 
centers of each segment. 
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Moore [12] examines this rediscretization technique as applied to the case of a uni- 
form circular vortex sheet in two-dimensions which has a known solution. For this flow, 
where there is no sheet stretching, the major concern is Kelvin-Helmholtz instability. 
Moore shows that the straightforward point vortex model experiences a discrete form of 
Kelvin-Helmholtz instability where the observed growth rate is that of the most unsta- 
ble mode. He then shows that the rediscretization technique of Fink and Soh acts as a 
smoothing mechanism which reduces the instability and delays the onset of the chaotic 
motion (but does not eliminate it). 

While these remarks on the error in the velocity are strictly applicable only to the 
two-dimensional sheet, it is assumed that the same general situation holds true for the 
three-dimensional axisymmetric sheets being considered here. The axisymmetric vortex 
segments behave locally as though they are two-dimensional, therefore, they should 
experience similar errors. To eliminate this error the rediscretization approach will be 
used for the three-dimensional axisymmetric sheet also. 

5.2.2 Application to the Discrete Axisymmetric Model 

Following Fink and Soh, after each time step a continuous model of the sheet is 
reconstructed by passing a smooth curve through all of the vortex elements in sequence 
along a cut through the sheet in a plane passing through the central axis. Simulta- 
neously, a functional representation of the circulation distribution along the sheet is 
obtained-. Given these, the total arc length is determined from the center to the outer 
edge of the sheet which when divided by the number of vortex elements yields a new 
sheet segment length. Vortex elements are then placed at the centers of each of the sheet 
segments and their strengths are determined from the sheet strengh distribution. If the 
element is a vortex filament then it is simply placed on the curve in the center of the 
segment. If vortex panels are being used then the panel center is positioned to coincide 
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with the point on the curve in the center of the segment and the orientation of the 
panel is set such that the panel is tangent to the curve at the panel center. The panel 
width is set equal to the arc length of the segment. Observe that this method does not 
yield a model of the sheet where the panel edges are joined into a continuous structure 
except in the limit as the panel width approaches zero. The sheet is not considered to 
be composed of panels but, rather, the flat panels are used as approximations to the 
smoothly curved sheet segments for purposes of velocity calculation. Since the panel 
velocity is calculated at its center, which lies on the sheet, this method yields velocities 
on the sheet itself. 

Once the sheet has formed a roll-up of about one and one half loops a different 
rediscretization technique is used in the spiral region. Rather than counting off equal 
length pieces of arc for the new segments, the segment lengths sure chosen such that the 
individual vortex elements will line up with each other on successive layers in the radial 
direction. This is done in an attempt to minimize the effect of the artificial singularities 
associated with each of the vortex elements. By arranging the elements in this manner, 
the effect of the singularities nearly cancels and the segment lengths of neighboring 
elements will still be very nearly equal. The remaining sheet segments in the flat part 
of the sheet are rediscretized with equal spacing as described before. 

The geometry and circulation of each portion of the sheet between any two vortex 
filaments or panel centers is described parametrically by an arc length parameter. For 
each ring interval t, the sheet is modelled as 


r% = r,(s) 

(5.30) 

Z{ = z,(s) 

(5.31) 

r< = r,(s) 

(5.32) 


where fj and Z{ represent the radial and axial components of geometry respectively, I\ 
is the circulation, and a is the arc length along the sheet along the cut starting from 
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the central axis. 


The functions for r,-, Z{, and T,- are a type of cubic spline that have the desirable 
property of bounding the interpolated curve which reduces the appearance of wiggles 
seen in other types of splines. The spline is constructed for a sheet segment between 
vortex elements t and t + 1 as shown in Figure 5.7. The first step is to construct a 
parabola, gi(s), through points * -2, * - 1, and * and another parabola, through 

the points t — 1, t, and i + 1. The interpolated function, /»•($), is then obtained as the 
following linear weighted average of gi(s) and 5<+i(s) 


/•(«) = <7i(s) 


g»+l ~ 3 

««'+! ~ Si 


+ 0t+i(s) 


a - Si 
*i + 1 ~ 


(5.33) 


where «,• and s,+i are the values of arc length at points t and $ + 1 respectively. This 
yields a cubic function for /,• which has a slope at the end points that matches that 
of the two generating parabolas. Repeated over all of the sheet segments, this process 
yields a continuous function from one end of the sheet to the other that has continuous 
first derivatives and is relatively free of the extraneous curvature seen in some other 
types of splines. 


The only exception to this scheme occurs at the outer edge of the sheet. As shown 
in Figure 5.8, the location of the outer edge of the sheet is not known. The known 
positions are at the vortex rings or panel centers, but the outer most ring or panel 
center is not placed at the outer edge of the sheet, it is at or near the middle of the 
final segment. For purposes of rediscretization the outer edge point of the sheet must 
be located. It is found by assuming that the final portion of the sheet has the shape 
of Kaden’s spiral, i.e., r = & 0 -2 / 3 , where r and 9 are referenced to the position of the 
center of the rolled-up tip vortex ring and & is chosen to match the Kaden spiral to 
the last known point on the sheet. As shown in Figure 5.8, the origin for 9 is taken 
as the horizontal axis through the spiral center and computed such that r(9) is not 
multivalued. 
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Sheet with 
inner spiral model 


Central vortex filament 
and higher order terms 



Figure 5.1: Inner Spiral Model 
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Figure 5.3: Numerical Integration over the Inner Spiral 



Figure 5.4: Rotation of Inner Spiral Model to Match Actual Sheet Orientation 
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Sheet at t = to 



is modelled as: 


tip vortex ring with finite core 
and higher order terms 
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Figure 5.5: Initial Representation of Sheet 
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Before 

Amalgamation 



After 

Amalgamation 


/ci, * 2 , and Kt are amalgamated to form k 


i 

t 


1. Strength of tip filament found by conserving circulation 

2. Location of tip filament found by conserving impulse 

3. Tip filament core radius found by conserving kinetic energy 


Figure 5.6: Sheet Edge Filament Amalgamation Procedure 


s = arc length 

fi = function being splined, either r,-, z i} or I\ 


/»■(*) = 9i(s) 


5 «+l ~ 3 
3 i+l ~ s i 


+ 9i+l{ 3 ) 


S — Si 
s «+l “ Si 


9i(s) = parabola through points t - 1, », » + 1 
= o 0 + ai(s - Si) + a 2 (s - s,) 2 

< 7 ,+i(s) = parabola through points i, t + 1, * + 2 
= b o + bi(s - s i+1 ) + b 2 [s - s i+1 ) 2 


Figure 5.7: Cubic Spline for Sheet Rediscretization 




Cubic Spline 



Figure 5.8: Rediscretization at Outer Sheet Edge 
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Chapter 6 

Application of Discrete Vortex Model 

6.1 Summary of Model 

1. An axisymmetric sheet is modelled as a series of equally spaced concentric vortex 
rings or panels. For vortex rings, the core diameter is chosen to be one half of 
the distance between consecutive filaments, this choice of core diameter conserves 
kinetic energy. 

2. The simulation starts at time t = to with the spiral portion of the tip replaced with 
the inner spiral model consisting of an equivalent single vortex filament and higher 
order terms. Strength, location, and core radius of the tip ring are determined from 
conservation of circulation, centroid of vorticity, and kinetic energy respectively. 
A typical non-dimensional value of to is 0.0075. 

3. The induced velocities are calculated using Equations 2.15, 2.16, and the inner 
spiral model of Section 5.1.2 and the positions of the elements are updated using 
a fourth order Runge-Kutta time step method. 

4. After each time step the sheet is rediscretized so that the vortex elements will 
remain equally spaced or along rays in the spiral region. 

5. Excess spiral loops are amalgamated into the tip ring vortex as they form. Strength, 
location, and core radius of the new tip ring are determined from conservation of 
circulation, centroid of vorticity, and kinetic energy. 
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6.2 Results 


The validity of the method is demonstrated by computing the roll-up of a two- 
dimensional sheet similar to that shed by an elliptically loaded wing. The results are 
compared to the similarity solution given by Kaden for the inner part of the roll-up 
of a semi-infinite sheet. This case is also used to compare and contrast the panel and 
filament models. Both are used to compute the two-dimensional roll-up and the results 
are compared. 

The panel model is then applied to vortex disks of two vorticity distributions, an 
elliptic distribution produced by a disk which translates through a fluid at constant 
speed and then suddenly dissolves, i.e., the three-dimensional analogue of the elliptic 
wing problem, and a distribution more typical of that for a helicopter rotor blade with 
the load increasing towards the tip. 

In all of the results that follow, t, T, length, and all velocities have been non- 
dimensionalized by i re /, To, b, and U rt f respectively, where, t re / = 6 2 /ro, b is the 
half-width of the sheet at time t = 0, To is defined in the loading equations that follow, 
and U re f is the magnitude of the velocity on the sheet at time t — 0. For comparison, it 
is recalled that Pullin’s non-dimensional time is given in terms of his similarity variable 
as r = A" 3 = a 4 i/r 3 . 

6.2.1 Two-Dimensional Model Validation 

The validity of the model may be evaluated by using it to compute the roll-up of a 
two-dimensional vortex sheet with an elliptic circulation distribution 

r = rWi-r* 
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which yields a value for the uniform downwash of U re f = Fo/26 and a corresponding 
sheet strength 

where r is the distance along the sheet. This is perhaps the most widely used test case 
for new methods because it is a good approximation to the three-dimensional roll-up of 
an elliptically loaded wing. Additionally, for very small values of time, an exact solution 
is available for the inner part of the spiral region in the form of a similarity solution. 
In the limit, as the center of the spiral is approached, the similarity solution becomes 
exact. Both the two-dimensional panel and filament models of Chapter 2 Me used and 
the results are compared. 

Simulation Parameters 

The parameters used in this simulation are as follows: 


Number of panels or filaments: 240 

Portion of sheet edge initially rolled: 2.5% 

Time step method: Runge-Kutta 

AXmax. 0.003 

Number of time steps: 3000 

Maximum number of loops in spiral: 10 


Initial Model of Sheet 

The initial sheet geometry and sheet strength distribution are shown in Figure 6.1. The 
same initial condition applies to both the panel and the filament models. Two and one 
half percent of the outer portion of the sheet rolled into the tip vortex filament results 
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in a tip fil am ent strength of about 22% of the total sheet strength. The large value is 
the result of the high concentration of vorticity in the outer portion of the sheet. As 
determined using the methods of Section 5.1.3, the location of the tip vortex filament 
is found to be r = 0.987 and z = 0.021 and the corresponding initial time is found to 
be to = 0.0123. As mentioned earlier, the vortex elements which represent the sheet 
initially lie flat horizontally and are stretched out as indicated in Figure 5.5. The initial 
slight curvature of the sheet is ignored. The corresponding sheet circulation is also 
linearly stretched along with the filaments. 

The initial vertical velocity distribution, v, is shown for the panel and filament 
moodels in Figure 6.2 and the initial horizontal velocity distribution, u, is shown for 
the two models in Figure 6.3. The exact initial condition for the vertical velocity on 
the sheet at time t = 0 is a uniform downwash from the center of the sheet to the outer 
edge and is equal to -1.0. The horizontal velocity should be zero everywhere except 
at the outer edge where both velocity components are undefined. The velocity profiles 
of Figures 6.2 and 6.3 are not expected to agree exactly with the true initial condition 
at time t — 0 since the initial condition for the numerical simulation is at a short time 
later, i.e., at t = to after a spiral has formed at the tip. However, to being small implies 
that the initial velocity profiles should be very similar to those at t = 0. For both 
panels and filaments the vertical velocity agrees very closely with the expected value of 
—1.0 over most of the sheet. The large values realized at the outer edge are attributed 
to the inability of either the panel or the filament models to accurately compute the 
induced velocity in a region where the vorticity gradients are very large, however, the 
panel model shows less error than the filament model. The initial horizontal velocity 
distributions for both panel and filament models are identical from the sheet center to 
the outer edge where the effect of the tip filament becomes apparent. Over most of 
the sheet, the horizontal velocity is zero but at the outer edge it increases sharply in 
the outward sense due to the tip vortex producing great stretching in this region of the 
sheet. 
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Roll-Up Before Amalgamation 


The resulting sheet roll-up computed using the inner spiral model of Section 5.1.2 is 
presented at several times up to the point when the first ten turns of the spiral are 
formed in Figure 6.4 for the panel model and in Figure 6.5 for the filament model. 
Note that in the figures some of the curves have been displaced vertically from their 
true positions in order to separate them for clarity. In general the models yield a very 
smooth roll-up with no sign of instability or other commonly observed chaotic behavior. 
The only noticeable difference between the panel and filament models at this point is 
that the panel model appears to roll-up slightly faster. Ten turns are achieved for the 
panel model at a time of t = 0.122 while for the filament model ten turns have not 
formed until t = 0.136. This difference is due to the retarding effect of the concentrated 
vortex at the control point on the neighboring turn. 

An enlargement of the spiral region at the time when ten turns have formed is given 
in Figure 6.6 for the panel model and in Figure 6.7 for the filament model. Aside from 
the slight difference in size due to the difference in time, there is very little observable 
difference between the models. The use of the equi-angular rediscretization scheme is 
apparent in the inner part of the spirals but note that the outer most turn is composed 
of segments of equal length rather than equal angles. It was found necessary during the 
course of the calculations to rediscretize the outer turn of the spiral with equal length 
segments because the equi-angular method, together with the rapidly increasing radial 
distance of the sheet from the spiral center, resulted in sheet segments that were too 
unevenly spaced on this final turn. Unequal spacing of the vortex elements has been 
established as one of the causes of the breakdown in these methods due to the associated 
large errors in the velocity calculations. Therefore the equi-angular rediscretizaton could 
only be used on the inner part of the spiral where the distance of the sheet from the 
spired center increases only slowly. 
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The sheet strength distributions plotted against arc length, s, measured from the 
sheet center are presented in Figures 6.8 for the panel model and in Figure 6.9 for 
the filament model. The stretching of the outer edge of the sheet is reflected in the 
sudden reduction in sheet strength as the edge is approached. In comparing panels and 
fil am ents, the results are very similar with, perhaps, the sheet strength distribution for 
the panel model slightly smoother near the end of the sheet. It is not clear if the wiggles 
are real or if they are due to numerical inaccuracy. Some waviness is expected in that 
the spired is immersed in a strain field imposed by the remainder of the sheet. 

The vertical velocity of the vortex filament which represents the spiral center, v tip , 
is given in Figure 6.10 for the panel model and in Figure 6.11 for the filament model 
along with the value obtained from Pullin’s results for the semi-infinite sheet. Pullin 
presents the position of the tip in terms of the similarity variable from which the tip 
vertical velocity may be obtained as 

= Jt [°- 489 ( a 0 2/3 ] - 10 (6.1) 

where the constant —1.0 is included to yield the velocity of the tip in the absolute 
reference frame used in the simulations and a is the sheet strength constant presented 
in Section 5.1.2. For the two-dimensional elliptically loaded sheet, a = Tq/ \/2. 

The oscillations between times t = to and t = 0.065 are not believed to be real but 
only a result of the inexact nature of the initial representation of the sheet. They are 
caused by the uneven distribution of vorticity around the tip filament as the sheet rolls 
up around it during the first few turns of the motion. The slightly faster roll-up of 
the panel model is reflected in the higher frequency of the oscillations in Figure 6.10 as 
compared to Figure 6.11. 

The inner spiral model introduced in Section 5.1.2 could conceivably remove all of 
the oscillation but apparently the inexact nature of the initial representation of the 
sheet is enough to introduce the errors realized here. The model acts to considerably 
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lessen the amplitude of the oscillations but does not completely eliminate them. The 
effect of the model on the oscillations may be evaluated by comparing the simulation 
with the model in place (i.e., Figure 6.10) to the same simulation without the model. 
The tip filament vertical velocity for the first few turns of the roll-up using panels is 
compared with and without the inner spiral model in Figure 6.12. The effect of the 
model is to reduce the amplitude of the oscillations by 50% or more and to increase the 
rate of roll-up, as is seen indirectly from the frequency of the oscillations. 

The distribution of circulation as a function of radius from the spiral center is eval- 
uated by summing the strengths of the vortex elements enclosed by circles of increasing 
radius centered at the tip vortex filament. This process is illustrated in Figure 6.13. As 
the radius increases, eventually all of the vortex elements of the spiral are included by 
the radius f. Further increase includes the unrolled part of the sheet as well. 

The radial distribution of circulation for the discrete models is presented along with 
Pullin’s result in Figure 6.14 for the panel model and in Figure 6.15 for the filament 
model. Consistent with the result for the semi-infinite sheet for small values of time, the 
circulation, once rolled up, remains at a constant radius from the spiral center. This is 
indicated by the clustering of lines for different values of time. Pullin’s results yield a 
value of circulation about 5% smaller than those obtained here at the smaller values of 
radius, i.e., both discrete models predict that the circulation distribution for the elliptic 
loading is slightly more concentrated in the inner part of the spiral than that for the 
semi-infinite sheet. At larger values of f the nearly parallel curves correspond to the 
part of the sheet which has not yet become rolled up. Since a point vortex is used to 
represent the center of the spiral, the curve of T(f) begins at a non-zero value and is 
constant until the turns of the spiral are encountered. 

The total circulation contained in the spiral region, T, p , (somewhat arbitrarily de- 
fined as where the sheet first becomes vertical) is presented in Figure 6.16 for the panel 


69 


model and in Figure 6.17 for the filament model along with Pullin’s results. The panel 
and filament models yield very similar results and both are very close to Pullin’s result. 
This is to be expected since for small values of time when only the outermost portion 
of the sheet is rolled up, the loading for the semi-infinite sheet and the finite sheet are 
similar. By the time ten turns have formed, approximately 42% of the total circulation 
is contained in the spiral 

The distance of the sheet from the spiral center plotted against angle, r(0), is given 
in Figure 6.18 for the panel model and in Figure 6.19 for the filament model at the 
time just before the first amalgamation. Also given is Kaden’s similarity solution using 
Pullin’s constants for the inner part of the spiral. The distance for the panel model is 
everywhere smaller than the values obtained from Pullin’s results while for the filament 
model the opposite is true. The small oscillations seen in both the panel and filament 
results indicate the slight ellipticity of the sheet due to the imposed strain field. 

Roll-Up After Amalgamation 

The roll-up continues after the time when ten turns have formed to about time t = 8.5 
with the amalgamation procedure being used to limit the spiral to ten turns. Ten turns 
was selected as the cut-off for the spiral after several trial runs indicated that errors 
would be introduced into the solution (detected by tracking the total fluid impulse) 
after this point if the turns were allowed to continue to accumulate. This is thought 
to be due to lack of resolution and the close spacing of the inner layers of the spiral. 
However, in some cases, when the roll-up was continued without any amalgamation and 
without regard for the errors, spirals of considerably more turns were obtained before 
the sheet crossed over on itself. In addition, the panel model would allow for, typically, 
one or two more turns than the filament model. 

The sheet geometry for the panel model is presented in Figures 6.20 and 6.22 and 
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for the filament model in Figure 6.21 and 6.23 for a series of times between t = 1.09 
and t = 8.5, i.e, between the first amalgamation and the end of the simulation. The 
corresponding sheet strength distributions are given in Figure 6.24 for the panel model 
and in Figure 6.25 for the filament model. The roll-up continues smoothly, still with no 
sign of Kelvin-Helmholtz instability. At time about time t = 6.8 the spiral center (not 
the centroid of vorticity) has moved its greatest distance inward to r = 0.778 which is 
very near the expected asymptotic value of r = 0.785. The effect of the amalgamation is 
clear from the figures as the inner parts of the spiral gradually disappear, being absorbed 
by the tip filament. The sheet strength distributions have more of the small bumps seen 
earlier indicating the presence of the strain field due to the rest of the sheet. At this 
later stage in the roll-up, there is very little difference in the results between the panel 
and the filament models although there are some slight differences in the sheet strength 
distributions. 

The vertical velocity of the tip filament has approached the correct asymptotic value 
for this flow of v = —0.205 at about time t = 1.3 as shown in Figure 6.26 for the panel 
model and in Figure 6.27 for the filament model. The total circulation contained in the 
spiral is shown in Figures 6.28 and 6.29 against time for the two cues. Also shown is 
Pullin’s result for the semi-infinite sheet. By the end of the simulation, about 96% of 
the circulation has become rolled up, as the spiral has been defined. 

The radial distribution of circulation in the spiral is presented in Figure 6.30 for 
panels and in Figure 6.31 for filaments along with Pullin’s result. For larger values of 
time, in contradiction to the results observed for small values of time, the circulation 
distribution is not completely constant. As the sheet rolls up, the vorticity in the spiral 
layers is slowly drawn in towards the center increasing the concentration of vorticity 
in the spiral core. This indicates that kinetic energy is not being conserved in the 
calculation and this was confirmed by numerical calculations. 
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A measure of the accuracy of the numerical simulation may be obtained by checking 
for the constancy of the total fluid impulse and kinetic energy. Since these are invarients, 
any deviation indicates numerical errors. For the panel model, the impulse increased 
by about 0.19% over the duration of the simulation or by about 0.000063% per time 
step. The corresponding values for the filament model are an increase of about 0.25% 
for the entire simulation and 0.000083% per time step indicating the greater accuracy 
of the panel model. A detailed analysis of the simulations revealed that about half 
of the increase in impulse is due to inaccuracy in time stepping and the other half is 
due to the rediscretization procedure. As was observed in the discussion on the radial 
distribution of circulation, kinetic energy was not conserved as well as impulse. As 
shown by Batchelor [2] , the part of the kinetic energy which depends on the positions 
of the vortices is given by 

T = «<«j log r.y 

i * 

where the summation is over all combinations of t and j except when * = j. Over the 
course of the simulation to t = 8.5 the kinetic energy decreased by about 6.0% for each 
of the discrete models, or by about 0.002% per time step. 

Validity of Model 

These results indicate that the model performs acceptably, both with panels and fila- 
ments, allowing for the roll-up of a good number of turns in the spiral and over a long 
period of time. The small change in impulse indicates that the time stepping and redis- 
cretization procedures perform reasonably well although there is some error, noticeably 
in the kinetic energy which should be corrected in future work. The agreement with 
Pullin’s results for the semi-infinite sheet is in most cases good with the exception of 
the initial motion of the spiral center. Perhaps a better initial definition of the unrolled 
part of the sheet would eliminate the oscillations noted in Figures 6.10 and 6.11 for 
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the vertical velocity of the spired center. The results show that the model can be used 
to discern the slight ellipticity in the sheet shape due to the imposed strain field and 
to give useful information on the rate of roll-up and circulation distribution for other 
vortex sheets undergoing roll-up. 

Comparison of Panels to Filaments 

The primary advantage of panels over filaments is to moderately increase the accuracy 
of the calculations and to slightly increase the length of time that the spiral can roll-up, 
i.e., the number of turns that can form around the core. In tried runs, typically one 
or two additional turns could form in the spiral with panels before the sheet would 
cross over on itself. These benefits are attributed to the more accurate velocity calcu- 
lation provided by the panel formulas, especially when the panels approach each other 
closely. However, because the distance between the spiral turns decreases rapidly as the 
center of the spired is approached, the improvement seen with panels is small in terms 
of the number of additional turns that can form. In addition, some of the potential 
improvement of panels over filaments is not realized because of the sheet segment line 
up process that is used in the rediscretization procedure (Section 5.2.2). The line up 
feature eliminates much of the difficulty with the singularities associated with the vortex 
filaments because of cancellation due to the symmetry which is introduced. Therefore 
the potential advantage of using panels is not fully realized. 

Axisymmetric Configurations 

The following calculations are performed on axisymmetric vortex sheets. The panel 
model is applied to vortex disks of two vorticity distributions, an elliptic distribution 
produced by a disk which translates through a fluid at constant speed and then sud- 
denly dissolves, i.e., the three-dimensional analogue of the elliptic wing problem, and a 
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distribution more typical of that for a helicopter rotor blade with the load increasing 
towards the tip. 

6.2.2 Elliptic Loading Case 

The vortex panel model is used to compute the roll-up of a circular disk of vorticity 
with an elliptic circulation distribution. This is a convenient test case for two reasons. 
First, it is the axisymmetric analogy of the frequently studied two-dimensional elliptic 
wake that was discussed in the previous section. Second, some features of the final 
rolled-up geometry sure known. As described in Section 5.1, the edge of the vortex sheet 
is modelled as a single vortex ring of uniform vorticity which will eventually acquire all 
of the vorticity in the sheet due to the amalgamation process. The final characteristics 
of the ring may be determined by equating the strength, impulse, and kinetic energy 
of the ring to those of the initial disk as shown by Taylor [22]. This, of course, is not 
intended to imply that a single vortex ring of uniform vorticity is the correct final state 
of the roll-up process. Taylor finds that for the disk of elliptic loading with radius 6, 
translating at a speed U re f, the corresponding vortex ring would have a radius of 0.8166 
and a translation velocity of 0.436tf re /. These values should provide a useful check 
towards the asymptotic limit of the simulation. 

Simulation Parameters 

The parameters used in this simulation are as follows: 

Number of panels: 

Portion of sheet edge initially rolled 
Time step method: 

Ax max : 


200 

2.5% 

Runge-Kutta 

0.003 
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Number of time steps: 

Maximum number of loops in spiral: 


1100 

10 


Initial Model of Sheet 


The initial configuration of the sheet is that which would be obtained from a circular 
disk translating normally to its axis which suddenly dissolves leaving a circular vortex 
sheet. The resulting circulation distribution is the same as that for the previous two- 
dimensional case 


r = r 0 Vi - r 2 


which yields a value for the uniform downwash of U Tt f = flTo/46. The vorticity distri- 
bution across the sheet, obtained from the circulation as 


r 


= r 0 


-r 

vi -r 2 


is illustrated in Figure 6.32 along with the initial positions of the vortex panel centers. 
Again, the two-dimensional tip model is used to define the initial conditions for the 
axisymmetric roll-up since for small values of time the axisymmetric sheet edge behaves 
similarly to the two-dimensional sheet. As for the two-dimensional cue, the initial 
conditions yield an equivalent initial time of to = 0.0123 and a tip filament containing 
22% of the total strength of the sheet. The location of the tip vortex filament is r = 0.987 
and z = 0.021 above the flat part of the sheet. The corresponding initial axial and radial 
velocity distributions, v and u, are given in Figure 6.33. 


Roll-Up Before Amalgamation 

The sheet roll-up computed by the vortex panel model prior to the point when 10 turns 
have formed in the spiral is presented in Figure 6.34 for times t = 0.039, 0.077, and 
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0.126. The sheet has formed 10 turns around the core at time t = 0.126 at which 
point the amalgamation process begins. As for the two-dimensional case, the sheet 
rolls up very smoothly exhibiting none of the typical chaotic behavior discussed earlier. 
A detailed view of the spiral region at time t = 0.126 is given in Figure 6.35. Com- 
parison with Figure 6.6 reveals very little qualitative difference in the roll-up between 
the two-dimensional and axisymmetric sheets at this eaxly stage. Note that the non- 
dimensionalization used allows for direct comparison between the two-dimensional and 
axisymmetric results in the non-dimensional time variable t. 

The corresponding sheet strength distributions are given in Figure 6.36. As for the 
geometry, a comparison with the two-dimensional results show no notable differences. 
It is interesting to note the sharply defined location of the maximum sheet strength. 
The sheet inboard of this point shows , very little stretching while outboard the sheet is 
stretched by a great amount due to the proximity of the spiral. The distance of this 
point from the edge of the sheet is found to vary approximately as y/t. 

The axial velocity of the tip vortex ring is found to experience the same slight oscil- 
lations during the time when the first few turns are forming as did the two-dimensional 
sheet, presumably for the same reasons. As can be seen in Figure 6.37 the oscillations 
subside after the first 3 or 4 loops but then reappear at time t = 0.095. The cause of 
the second set of oscillations is unknown, but as will be seen, they die out quickly. 

The distribution of circulation within the spiral region during this period is presented 
in Figure 6.38 along with Pullin’s result for the semi-infinite sheet. The constancy of 
the distribution with time is less for the axisymmetric case than it was for the two- 
dimensional sheet. At least part of this is due to the slightly irregular motion of the 
tip vortex filament since the distance plotted in the figure is measured from the tip. As 
for the two-dimensional case, Pullin’s results indicate a slightly less concentrated spiral 
core than that obtained here for the axisymmetric sheet. 
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The total circulation in the spiral region as a function of time is presented in Figure 
6.39 along with Pullin’s results. The spiral region is defined to start where the sheet 
first becomes vertical. For the initial part of the roll-up, the agreement is very close 
between Pullin’s results and the axisymmetric discrete model indicating that the initial 
roll-up is very nearly two-dimensional. 

The distance of the sheet from the spiral center, f, as a function of the angular 
position 0 at time t = 0.126 is given in Figure 6.40 along with Pullin’s results applied 
to Kaden’s spiral. The small oscillations of the numerical results about the Pullin value 
indicate a slight ellipticity in the sheet shape. This is assumed to be caused by the strain 
field imposed on the spired from the rest of the sheet. The sheet strength distribution 
as a function of 9 is given in Figure 6.41 at time t = 0.126. 


Roll-Up After Amalgamation 

The roll-up calculations are continued past time t = 0.126 by using the amalgamation 
procedure to limit the number of turns in the spiral to 10. As for the two-dimensional 
case, 10 turns was selected as a compromise between achieving good resolution in the 
spired and suffering from numerical inaccuracy due to the vortex elements approaching 
each other too closely in the inner turns. The sheet geometry and strength distributions 
are presented in Figures 6.42 and 6.43 respectively for times t = 0.32, 0.65, and 1.52 
when the calculation was stopped. The results are characterized by smooth roll-ups with 
no chaotic behavior. A feature of the axisymmetric roll-up that is not present in the 
two-dimensional case is the development of the asymmetry in the spiral which becomes 
apparent by time t = 1.52. The spiral has taken on a distinctly lopsided appearance 
with the outer most turn expanding inward and upward. Examination of the sheet 
strength distributions reveals the development of a local concentration of circulation in 
that part of the sheet. This is undoubtedly a function of the axisymmetric nature of 
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the sheet since it is not present at all in the two-dimensional case. 


The axial velocity of the tip ring, t > t , p , is shown in Figure 6.44 along with the 
asymptotic value given by Taylor of v = -0.436. After the previously mentioned second 
set of oscillations dies out, the asymptotic value is approximately reached at time t = 
1.02. At this point, about 73% of the circulation has been rolled into the spiral as can 
be seen in Figure 6.45, where the spiral circulation is plotted against time. By the end 
of the simulation at time t = 1.52, about 80% of the circulation has been absorbed by 
the spiral. 

The circulation distribution in the spiral as a function of radius is presented in 
Figure 6.46. The lack of symmetry in the spiral at large values of time is responsible for 
the jagged appearance of the curves. However, the inner turns of the spiral are more 
uniform and show a similar concentration of vorticity as time increases as was seen for 
the two-dimensional case. 

The distance of the sheet from the spiral center, f, and the sheet strength distri- 
bution, 7 , as functions of 9 are presented in Figures 6.47 and 6.48 respectively. The 
oscillations in the plot of f(9) are much larger compared to the earlier time indicating 
that the spiral turns are becoming more elliptic, especially the outermost turns. 

The fluid impulse is seen to vary by about 1 % over the duration of the simulation, 
or, about 0.001% per time step. This is about 5 times larger than the error realized in 
the two-dimensional case and is attributed to the approximations made in the derivation 
of the velocity formulas for the vortex panels. 
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6.2.3 Helicopter Loading Case 


Configuration Description 

A vortex disk similar qualitatively to that which would be produced by a single heli- 
copter rotor blade with the load increasing towards the tip is modelled as an example of 
a more complex vorticity distribution. The rotor blade loading and resulting vorticity 
distribution across the sheet are given by 

T(r) = IWTT7* (6.2) 

-r(r) = r 0 £ (r (6.3) 

Note that inboard of r = 0.82, the sheet strength is positive and outboard of r = 0.82 it 
is negative with the total circulation of the sheet equal to zero. In a gross sense, this is 
equivalent to two coplanar vortex rings of equal but opposite strength at different radii. 
As will be seen, the resulting sheet motion agrees with this interpretation. 

The nature of the roll-up for this cue is fundamentally different from that of the 
elliptic loading case in that there are two regions on the sheet of high vorticity concen- 
tration. The first is at the edge of the sheet as for the elliptic loading and the second 
is inboard at a radius of about r = 0.52. Apparently these differences in sheet strength 
distributions have resulted in some undesirable effects in the roll-up calculations. The 
sheet repeatedly experienced instabilities at the inner point of maximum sheet strength. 
These were characterized by an oscillation in the sheet strength distribution near the 
point of maximum strength with a period of about 9 or 10 panel widths. A similar 
phenominon was observed by Sugioka [21] where the period of oscillation was about 
11 panel widths. It was found that the instability could be delayed significantly by 
replacing the cubic fit of the circulation in the rediscretization procedure with a linear 
function. In addition, increasing the number of vortex panels tends to aggravate the 
problem. For example, using 200 panels to model the sheet results in the sheet becoming 
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unstable before the secondary inner vortex can begin to form. Since this is one of the 
more interesting features of the flow, it was necessary to reduce the number of panels 
from that for the elliptic loading case to increase the length of time that a valid roll-up 
could be observed. For this reason, the calculation was limited to 160 panels. Also, 
because of the much greater lengthening of the sheet due to the stretching realized in 
this case, only four spired loops are allowed to form around the tip ring in order to 
minimize the total sheet length and maximize the panel density along the sheet. 

Simulation Parameters 

The parameters used in this simulation are as follows: 


Number of panels: 160 

Portion of sheet edge initially rolled: 2.5% 

Time step method: Runge-Kutta 

Axnu,,: 0.003 

Number of time steps: 1900 

Maximum number of loops in spiral: 4 


Initial Model of Sheet 

The initial sheet geometry and sheet strength distributions sure shown in Figure 6.49. 
Two and one half percent of the outer portion of the sheet accounts for about 55% of the 
negative vorticity leaving the sheet itself to be of mostly positive vorticity. The location 
of the tip ring is determined as for the previous case, by applying Kaden’s solution to the 
outer part of the sheet, although its validity is somewhat more questionable here since 
the helicopter circulation distribution is even further removed from the parabolic one of 
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the semi-infinite sheet than that of the elliptic loading. The radial and axial location of 
the tip vortex ring are the same as for the elliptic case, r = 0.987 and z = 0.021 above 
the plane of the flat part of the sheet which is taken as z = 0 at t = to- The initial time 
is also the same, to = 0.0123. 

The initial axial and radial velocity distributions, v and u, are shown in Figure 
6.50. Although a simple analytic expression for the initial downwash distribution is not 
available in this case, the computed results at least appear as would be expected with 
the interpretation of the sheet mentioned before as two coplan ar vortex rings. The inner 
part of the sheet is moving upward and the outer put is moving downward. The initial 
radial velocity of the sheet is very similar to that of the previous cue as is expected 
since the properties of the tip ring for the two cues are nearly the same. 

Sheet Roll-Up 

The sheet geometry for a series of times to t = 3.51 is presented in Figures 6.51 and 
6.52 and the corresponding sheet strength distribution in Figures 6.53 and 6.54. The 
resulting sheet motion yields two distinct regions on the sheet. The outer region shows 
a roll-up very similar to that for the elliptically loaded disk especially up to a time of 
t = 0.66. This is not unexpected in that the sheet strength is very similar for the two 
cues near the edge. Four spiral turns, the maximum, are formed at time t = 0.076 after 
which the amalgamation procedure lumps any additional turns into the tip filament. 
The inner region of the sheet initially shows a smooth upward motion opposite that of 
the outer part of the sheet. This is due to the oppositely signed vorticity in this region. 

As the roll-up progresses, the part of the sheet with the positively signed vorticity 
begins to move into the outer layers of the spiral. This is made evident by the swelling 
of the outer turn starting at time t = 0.66. By the end of the simulation the oppositely 
signed vorticity hu caused the outer turn to move away from the spiral considerably, 


81 


especially in the inboard upper quarter. The amalgamation process resulted in the tip 
vortex filament absorbing about 95% of the negative vorticity in the sheet. 

The most interesting feature of the flow is the development of the secondary con- 
centrated inner vortex structure seen to be forming at a radius of about r = 0.45 and a 
height of z = —0.15 at time t = 1.94. Comparison with the sheet strength distributions 
indicates that the location of the inner vortex corresponds to the very strong maximum 
observed at an arc length of a = 0.95. As time increases the maximum becomes more 
and more sharply defined until it consists of only a few of the vortex panels. This 
maximum is a remnant of the original peak in the positive vorticity distribution at time 
t = 0 and occurs at very nearly the same material particle. The position of the material 
particle which was at the peak at time t — 0 is indicated in the sheet strength figures by 
the large solid dot. By the end of the simulation, this particle has moved away from the 
peak by only two or three panel widths. Although it seems clear that the inner vortex 
structure would assume the configuration of a double branched spiral if the simulation 
could be run for long enough, it would take a very long time, i.e., the inner and outer 
roll-up procede at very different rates. By the end of the simulation, the inner vortex 
structure has formed about one half of a spiral turn. 

The calculations invariably break down shortly after this point as the sheet crosses 
over itself in the region of the inboard secondary vortex. This is assumed to be caused 
by the very high strength of a few of the panels as is seen in Figure 6.54 at time 
t = 3.51. The very high concentration of vorticity that is forming in the secondary 
vortex, spiral would perhaps be best handled by another concentrated filament as was 
done by Hoeijmakers and Vaastra [9]. 

Examination of the sheet strength at times t = 2.57 and t = 3.51 in Figure 6.54 
reveals the development of a potential instability in the sheet at a point just inboard 
of the sharp maximum coinciding with the inner vortex. The fluctuations of strength, 
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which become larger with time, are characteristic of the instability that has resulted in 
unsucessfitl runs with larger numbers of panels. However, in this case the disturbance 
grew slowly enough so that the simulation was able to continue to the point when the 
secondary inboard vortex could develop. 

The axial velocity of the tip filament, v t i p , is given in Figure 6.55. The inset in 
the Figure details the motion for early times. Initially the motion is very similar to 
that for the elliptic loading case but for later times the increase in speed is less due to 
the lesser total negative circulation in the spiral. This case shows more of the erratic 
behavior of the tip filament as the first few turns are formed than was seen in the 
elliptic loading case. This may be partly because the inner spiral model is based on a 
parabolic circulation distribution which is not very accurate away from the sheet edge 
in the helicopter loading case. Therefore for larger values of time the entrainment of 
the increasingly non-parabolic circulation distribution into the spiral is likely to have 
an effect on the inner spiral characteristics which is not accounted for in this model. 

The total circulation of the edge roll-up is shown in Figure 6.56 against time. The 
different nature of the loading for the helicopter case is apparent when comparing this 
Figure against the corresponding one for the elliptic loading case. Initially, the results 
are consistent with Pullin’s for the semi-infinite sheet as they should be considering 
that the sheet edges are similar. However, at time t = 0.83 the circulation in the spiral 
reaches a maximum at which point it contains all of the negative vorticity in the sheet 
and then begins rolling up the positive vorticity resulting in the total spiral strength 
being reduced. The spiral strength is reduced by about 12% from time t = 0.83 to 
t = 3.51 at the end of the simulation due to the addition of the positive vorticity. 
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Figure 6.4: Sheet Roll-Up Before Amalgamation - Panel Model 
Two-Dimensional Sheet 
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Figure 6.5: Sheet Roll-Up Before Amalgamation - Filament Model 
Two-Dimensional Sheet 


87 







Figure 6.8: Sheet Strength Distribution to i = 0.122 - Panel Model 
Two-Dimensional Sheet 



Figure 6.9: Sheet Strength Distribution to f = 0.136 - Filament Model 
Two-Dimensional Sheet 
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Figure 6.10: Vertical Velocity of the Spiral Center to t = 0.122 - Panel Model 

Two-Dimensional Sheet 
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Figure 6.11: Vertical Velocity of the Spiral Center to t = 0.136 - Filament Model 
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Figure 6.12: Effect of Inner Spiral Model on Tip Filament Vertical Velocity 




Figure 6.13: Definition of Radial Circulation Distribution 
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Figure 6.14: Radial Distribution of Circulation to t = 0.122 - Panel Model 

Two-Dimensional Sheet 
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Figure 6.15: Radial Distribution of Circulation to t = 0.136 - Filament Model 

Two-Dimensional Sheet 
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Figure 6.16: Circulation Contained in Spiral to t = 0.122 - Panel Model 

Two-Dimensional Sheet 
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Figure 6.17: Circulation Contained in Spiral to t — 0.136 — Filament Model 

Two-Dimensional Sheet 
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Figure 6.33: Axial and Radial Velocity Distributions at Time t = to - Elliptic Loading 

Axisymmetric Sheet 
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Figure 6.34: Sheet Roll-Up Before Amalgamation - Elliptic Loading 

Axisymmetric Sheet 


Figure 6.35: Detailed View of the Spiral Region with Ten Turns - Elliptic Loading 

Axisymmetric Sheet 
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Figure 6.36: Sheet Strength Distribution to t = 0.126 - Elliptic Loading 
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Figure 6.37: Axial Velocity of the Spiral Center Filament to t = 0.126 - Elliptic Loading 

Axisymmetric Sheet 
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Figure 6.38: Radial Distribution of Circulation to t = 0.126 - Elliptic Loading 

Axisymmetric Sheet 



t 

Figure 6.39: Circulation Contained in the Spiral to t = 0.126 - Elliptic Loading 

Axisymmetric Sheet 
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Figure 6.51: Sheet Geometry to t = 0.076 - Helicopter Loading 
Axisymmetric Sheet 
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Figure 6.52: Sheet Geometry to t = 3.51 - Helicopter Loading 
Axisymmetric Sheet 





Figure 6.54: Sheet Strength Distribution to t = 3.51 - Helicopter Loading 

Axisymmetric Sheet 
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Chapter 7 

Conclusions 


Vortex sheet motion is an interesting and challenging subject. As am initial value 
problem it is ill-posed as the wavelength of disturbances approaches zero and is generally 
plagued with instabilities. However, it is possible through smoothing techniques to 
calculate the motion and determine useful information at moderate scales about the 
subsequent roll-up. This work was directed at developing a reliable and accurate method 
for calculating the motion of axisymmetric vortex sheets, especially the roll-up at the 
sheet edge. Knowledge of this type of motion is useful in the study and design of aircraft 
of various types including those of fixed and rotary wings. 

The motion of the vortex sheet was studied in a Lagrangian manner, tracking in- 
dividual material particles on the sheet. Both vortex panels and vortex filaments were 
used to represent the sheet. The material particles at which the velocities are calculated 
are located at the panel centers or at the filaments. A complete set of equations for the 
induced velocity was presented for both panel and filament models. 

A comparison between the Euler and the Runge-Kutta time stepping procedures 
yielded the result that the Runge-Kutta procedure is far more accurate than the Euler 
procedure for equivalent step sizes in vortex motion calculations. Even considering the 
additional computations that must be performed at each time step for the Runge-Kutta 
method it still outperforms the Euler method by many times in vortex calculations. 

The inherent difficulties in the calculation of vortex sheet motion were revealed 
in a straightforward calculation of the motion of an elliptically loaded disk of vorticity. 
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Shortly after starting, the vortex filaments used to model the sheet exhibited the classical 
chaotic behavior seen by many other investigators. It was determined that inaccurate 
time stepping was not the cause but that a variety of problems associated with the sheet 
itself and the sheet discretization procedure were responsible. 

The largest single cause of the chaotic behavior observed very early in the simulation 
was determined to be due to the inability of a finite number of discrete vortex elements 
to model the singularity at the sheet edge. This was resolved with the introduction of 
a special model of the edge of the sheet. The model replaces the outermost portion of 
the sheet with a single vortex filament of equivalent strength and a number of higher 
order terms which account for the asymmetry of the spiral. 

Another difficulty associated with representing the vortex sheet with a series of 
discrete vortex elements is the error associated with non-uniform element separation. 
As the sheet stretches unevenly, the initially equally spaced vortex elements become 
unevenly spaced which can introduce severe errors in the velocity calculations. The 
solution to this problem wets to rediscretize the sheet at every time step to ensure that 
the vortex elements remain equally spaced. 

The validity of the model was tested by computing the roll-up of a two-dimensional 
elliptically loaded sheet. A smooth roll-up was observed for the entire period of the 
simulation indicating that the method appears to have eliminated the primary causes of 
the chaotic motion observed earlier. Comparing the results to those obtained by Pullin 
for the similar semi-infinite sheet demonstrate that the model may be used successfully 
to simulate the roll-up of a vortex sheet, although the motion for very early times shows 
some error in the motion of the spiral center. Additionally, the two-dimensional test 
case was used to compare the panel and filament models. Few differences were observed 
in the results between the panel and filament models. The primary advantage of panels 
over filaments is their increased accuracy resulting in the ability to form more turns in 
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the spiral. However, the line up of the vortex elements in the spiral served to eliminate 
much of the problem of the singularity associated with the vortex filaments and so the 
potential advantage of panels was not fully realized. 

The model was used to calculate the roll-up of two axisymmetric disks, one with an 
elliptic loading and one with a. loading similar to that of a helicopter rotor in hover. 
For the elliptically loaded disk the method yielded a very smooth roll-up with none of 
the instability or chaotic behavior seen in the straightforward case. Fluid impulse was 
conserved to within about 1 % indicating that the method is accurate. The helicopter 
loading case initially suffered some instability in the region of the inboard sheet strength 
maximum, but this was largely overcome by replacing the cubic fit of circulation with 
a simpler linear fit. The cubic function had introduced some extraneous curvature into 
the circulation distribution which became magnified by the differencing process used to 
obtain the strengths of the individual panels. The simulation of the helicopter rotor 
wake successfully captured the development of the secondary inboard vortex although 
it was not capable of fully simulating the double branch roll-up past the early stages. 
However, the addition of another isolated vortex at the inboard vortex position would 
likely allow the roll-up to continue. 
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Appendix A 

Vortex Polygon Model 


A.l Cut-Off Method 

The self-induced velocity of a vortex filament is infinite at any point where the 
curvature is non-zero. A popular method to avoid this problem when attempting to 
calculate the velocity on the curved filament by means of Biot-Savart integration, is to 
remove a section of the filament around the point of calculation from the integration. 
This is commonly known as the cut-off method. As shown in Figure A.l, a segment of 
length l e is removed from either side of the control point. Bliss [4] has shown that the 
correct self-induced velocity is obtained if the cut-off length satisfies 

log- = -log2+ ^ -A + C (A.l) 

CL 

where c is the radius of the equivalent core section of uniform vorticity. A and C are 
constants which are determined by the axial velocity and the distribution of vorticity 
within the core respectively. For no axial velocity, A = 0, and for uniform vorticity, 
C =1/4 which yields l e /c = 0.6420. 


A.2 Polygon Model 

The approach considered in this appendix is to model a vortex ring as a vortex 
polygon consisting of N straight sides of equal length. A vortex polygon is constructed 
from a circular vortex ring as illustrated in Figure A.2. The ends of the straight segments 
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lie on the ring and the segments themselves lie inside the circular ring. The motion of 
the polygon may be tracked by computing the velocity at the centers of each of the 
segments due to all of the other segments and all of the other polygons. This naturally 
results in a type of cut-off procedure where the cut-off length, / e , is half of the length of 
one of the polygon sides. The objective is to determine the appropriate cut-off length 
in order to model a vortex ring of a given core section radius. 

For a vortex ring of radius R, and strength k, the self-induced axial velocity is given 
by Lamb [11] as 



where e = c/R, the ratio of core radius to ring radius. Solving for e yields 

e = 8 exp|^- + V— J (A.3) 

The idea is to numerically evaluate the velocity at the center of one of the straight 
segments by applying the Biot-Savart integral to all of the other segments. The resulting 
velocity is then substituted into equation A.3 to determine an effective e for this cut-off 
length. 

For large values of N (i.e., N > 75), the value obtained is l c /c = 1.387. This is 
considerably larger than the cut-off length determined by Bliss because of the difference 
in the location of the control point. For Bliss the control point is located on the ring 
itself, while for the polygon method, the control point is located on the center of one of 
the segments which is further in towards the center. 

The implication of this result is that the modelling of a vortex line by a set of straight 
filament sections in which the self-induced velocity of the straight filament section upon 
itself is zero is equivalent to a vortex core of radius c = 1.387/ c . 

This result yields a method for determining the number of polygon segments to be 
used in a ring system to model a disk. Consider a circular vortex disk of radius 6 which 
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is subdivided into a number of equally spaced annuluses. At the center of each annulus 
a vortex polygon is placed having the same strength as the annulus which it represents. 
An effective core radius may be determined for the polygon using a result presented in 
Section 2.3, i.e., c = w/4 where w is the width of the annulus. If the cut-off length 
is indicated by l e and the associated subtended angle by 9 then for large values of N, 
i.e., small values of 9, l c may be approximated by l e /R « ir/N or with e = c/R and 
l e /c = 1.387 we get 

N = 1.766 Nr*- (A.4) 

0 

where Nr is the number of polygons used to model the disk. This expression yields 
the number of segments to use in each of the polygons from the center ( r/b — 0) to 
the outer edge ( r/b = 1). The number of polygon segments increases linearly with the 
distance from the center of the disk with the value at the center equal to zero. Clearly 
this would have to be approximated, but that is acceptable in that the strength of the 
rings in the center is very small so the error should also be small. Also, the number of 
polygon segments increases linearly with the number of polygons. This could represent 
difficulties computationally in that a very large number of polygon sides may be needed 
to model a vorticity disk and achieve good resolution. 
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Appendix B 

Computer Program ROLLUP 


B.l Description of Program 


All of the results presented in Chapter 6 were calculated by the computer program 
ROLLUP. ROLLUP incorporates the vortex dynamics algorithm presented in Chap- 
ters 2 through 5 including the velocity calculations, the inner spired model, and the 
rediscretization procedure. The operation of the program is detailed in the series of 
flowcharts presented in Figures B.l through B.5. 

The user supplied input data to program ROLLUP is described in what follows and 


is read by 

ROLLUP from unit 5. 

IWR 

= 

FORTRAN unit number for printed output. 

NRING 


Number of vortex elements to be used initially, the actual number 
used may vary according to the value supplied for MINSECT. 

IAXI 

= 

0 for two-dimensional sheet. 


= 

1 for axisymmetric sheet. 

INTKOD 

= 

0 for Euler time stepping procedure. 


= 

1 for Runge-Kutta time stepping procedure. 

DXMAX 


Ax max from Section 3.1, the maximum distance moved 
by any vortex element within one time step. 

NSTEP 

= 

Number of time steps. 

THMAX 

= 

The maximum angle that the sheet may roll-up. Any roll-up 
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FRAC 

IWING 

IREDIS 

NU2 

IPANEL 

LINEUP 

MINSECT = 

MAXSEGS = 


greater than THMAX is amalgamated into the tip filament. 
Fraction of the sheet initially rolled (see Section 5.1.3). 

0 for elliptic loading. 

1 for helicopter loading. 

Flag for allowing calls to subroutine REDISC. 

0 for rediscretization not performed. 

1 for rediscretization performed at every time step. 

Frequency of unformatted WRITEs to unit 2 for plotting. 

0 for filaments. 

1 for panels. 

0 to inhibit the spiral lineup option. 

1 to use the spiral lineup option. 

The minimum number of sheet segments to be used on each 
turn of the spiral during rediscretization. The total number of 
vortex elements may be increased during execution, if necessary, 
to satisfy MINSECT. This parameter has no effect if LINEUP = 0. 
The maximum number of vortex elements that the program 
is allowed to use when trying to satisfy MINSECT. This 
parameter also has no effect if LINEUP = 0. 


The main program, Figure B.l, acts primarily as a controller, calling the principle 
subroutines inside of the time step loop. First the user supplied data is read from unit 
5. This is followed by the program initialization. All variables are initialized including 
the vortex element positions and strengths and if panels are used, their orientations. 
The properties of the tip vortex filament are defined and the equivalent initial time 
is calculated. Next, all of the initial data, including the initial sheet velocities are 
written to unit 2. At this point the time step loop is entered where the three primary 


128 


subroutines, MOVEM, REDISC, and LUMPS are called. At the end of each time step, 
depending on the NU2 parameter, the resulting data is written to unit 2. 

Subroutine MOVEM updates the positions of the vortex elements on the sheet and 
the tip vortex filament. As illustrated in Figure B.2, the first step is to decide be- 
tween the Euler and the Runge-Kutta time step procedures based on the user input. 
If the Euler time step is used then Subroutine VELOCY is called which calculates the 
total induced velocity at every vortex element. This is followed by a call to subrou- 
tine UPDATE which moves the vortex elements according to Equations 3.1 and 3.2. If 
the Runge-Kutta time step has been selected then a similar process is followed except 
that the Runge-Kutta coefficients of Equations 3.3 and 3.4 are calculated by subroutine 
GETK after the velocities are determined by subroutine VELOCY. Then the vortex po- 
sitions are updated according to the Runge-Kutta formulas and the process is repeated 
using the modified vortex positions. The cycle is completed four times to generate the 
four sets of Runge-Kutta coefficients. The final positions are then obtained using the 
four Runge-Kutta coefficients in Equations 3.3 and 3.4. 

Subroutine VELOCY, outlined in Figure B.3, acts as a switch to select from the four 
possible velocity calculation subroutines: axisymmetric panels in subroutine PANVEL, 
axisymmetric filaments in subroutine INDVEL, two-dimensional panels in subroutine 
PANVEL2D, and two-dimensional filaments in subroutine INDVEL2D. 

Subroutine REDISC is responsible for rediscretizing the sheet at every time step. As 
illustrated in Figure B.4, the first step is to fit the cubic splines of Equations 5.30 and 
5.31 to the current vortex positions to reconstruct the sheet. From this, the total sheet 
length is calculated. If the lineup option is not being used, a new, constant, panel width 
for the entire sheet is found by dividing the total sheet length by the number of vortex 
elements. Then the sheet is rediscretized by counting off equal length segments along 
the sheet, each of width equal to the new panel width. If the lineup option is being used 
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then first a test is made to determine if there is at least one and one half turns in the 
spiral. If not, then the sheet has not rolled up sufficiently to use the lineup option and 
the sheet is rediscretized with the equal panel width method. If the spiral is sufficiently 
formed then the lineup process continues with the determination of the length of the 
sheet forming the turns of the spiral and the number of segments that can be put on 
each turn. If the MINSECT parameter is not satisfied then the program attempts to 
increase the total number of vortex elements in the sheet until it is. This, however, is 
subject to the MAXSEGS parameter which is an upper limit to the number of elements 
that may be used. The new vortex element positions are then found in the spiral such 
that the elements line up along rays emanating from the origin. The remaining part of 
the sheet is rediscretized with equally spaced elements. 

At this point, the sheet geometry has been rediscretized and all that remains is to 
determine the strengths of the vortex elements. The vorticity distribution along the 
sheet is integrated from the center of the sheet to the edge to yield the circulation 
distribution. Equation 5.32 is used to spline fit the circulation which is then differenced 
to obtain the individual vortex element strengths. 

Subroutine LUMPS performs the amalgamation procedure as illustrated in Figure 
B.5. First, the extent of the roll-up is calculated by measuring the angle of the final sheet 
panel to determine if amalgamation is necessary, if not, the subroutine is exited. If the 
sheet has formed enough turns then any element with an angle greater than THMAX 
is considered to be removed from the sheet and amalgamated into the tip filament. 
The circulation, impulse, and kinetic energy of the original sheet are found. The new 
tip filament strength is found from conservation of circulation. The position of the tip 
filament is found from the conservation of impulse and the new core section radius of 
the tip filament is found from the conservation of kinetic energy. Before LUMPS is 
exited, Subroutine REDISC is called to rediscretize the sheet and restore the number 
of vortex elements in the sheet to the value before the amalgamation. 
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Figure B.l: Main Program of ROLLUP 
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Figure B.2: Subroutine MOVEM 
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Figure B.3: Subroutine VELOCY 
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Figure B.4: Subroutine REDISC 
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Figure B.5: Subroutine LUMPS 
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B.2 Program Listing 


PROGRAM ROLLUP 

C — This program calculates the trajectories oi a group of vortex rings 
C or panels which simulate a vorticity disk. The user has the option of 
C using Euler or Runge-Kutta time integration. The curve fitting method 
C allows for the inclusion of a concentrated tip core which is not part 
C of the sheet. At any time step if the sheet rolls up by more than 
C a certain amount (THMAX) the excess part of the sheet is 
C chopped off and its vortcity lumped into the tip core. In 
C doing this, conservation of impulse determines the location 
C of the resulting core and conservation of kinetic energy 
C determines the core radius. 

C This version allows the simulation to start at some time t>0 by 
C lumping a fraction of the outer portion of the sheet into the 
C concentrated tip vortex at the start. 

C 

COMMON /Bl/ RADIUS (0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0 : 900). CORE (0:900), TIME (0: 5000), TENERG (0:5000), 

$ TIMPUL (0:6000) 

COMMON /B2/ NRING , PI .DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
C COMMON /B3/ DMIN (0:900) 

COMMON /TH/ THETA (0:900), THMAX. DT (0:900) 

COMMON /B9/ IWR 

COMMON /PANEL/ IPANEL , DRPAN (0 : 900) , DZPAN (0 : 900) . WPAN (0:900) 

COMMON /Cl/ LINEUP. MAXSEGS.MINSECT 

COMMON /D1 / DISTMX 

COMMON /SEGS/ AXLSEG (0 : 900) 

COMMON /FI/ ALPLAST , GMTO 

COMMON /ZZ/ ZXO.ZYO 

COMMON /Gl/ TIPVAX 

COMMON /E2/ DUM1 , DUM2 , DUM3 , DUM4 , DUM5 , AE2 

COMMON /NOW/ TIMENOW , B IGRREF 
COMMON /El / TAUO , BDUM , CDUM , DDUM , TAU2L00P 
DIMENSION SGAM(0:900) ,DDGM(0:900) 

REAL*4 LAMBDA. K 
CHARACTER PLTITL*80 , TITLE*80 
C 

C NRING ■ total number of rings (counting starts at 0) 

C NSMALL ■ number of vortex rings which form the sheet. NSMALL=NRING- 1 
C 

C Ring Indie ies 

C 

C Index Meaning 

C 

C 0 Degenerate ring on the central axis of the ring system. 
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C This ring has zero circulation and is used only to 

C define the inner boundary of the sheet. 

C 

C 1 - NSMALL These are the "real" vortex rings, i.e., they have non- 

C (inclusive) zero circulation. These all lie on the sheet. 

C 

C NSMALL This ring forms the outer boundary of the sheet. 

C 

C NRING This is the index of the very last ring. It is the 

C concentrated core tip ring. It does not lie on the sheet. 

C Initially, this ring has circulation determined by FRAC. 

C As the sheet rolls up, the circulation from the outer 

C portions, which is chopped off, gets put into the core. 

C 

PI-ACOS(-l.O) 

TIME(1)=0.0 

C 

WRITE (6,*) * ENTER UNIT NUMBER FOR PRINTED OUTPUT’ 

READ(5,*)IWR 

INEW-0 

WRITE(6, *) * ENTER: NRING, IAXI (1=AXI 0-2-D) * 

READ (5,*) NRING, IAXI 
ILOOl 

WRITE(6,*)’ ENTER: (1) FOR RUNGE-KUTTA OR (0) FOR EULER’ 

READ ( 5 , * ) INTKOD 

WRITE (6,*) ’ ENTER: DXMAX, NSTEP (IF DXMAX .LT. 0, DT=DXMAX) ' 

READ (5 , *) DXMAX , NSTEP 

WRITE (6,*)’ ENTER: THMAX (MAX ROLLUP ANGLE BEFORE CUT-OFF)’ 
READ(5,*)THMAX 

WRITE (6,*)’ ENTER: FRAC (OUTER PORTION TO BE LUMPED AT START)’ 
READ(6,*)FRAC 

946 WRITE(6,*) ’ ENTER: IWING (EXPONENT ON LOADING CURVE 0,2,4)’ 
READ(5,*)IWING 

IF (IWING. NE.O. AND. IWING. NE. 2. AND. IWING. NE. 4) GO TO 946 
WRITE(6,*) ’ ENTER (1) FOR REDISCRETIZATION (0) FOR NONE’ 
READ(5,*)IREDIS 
WRITE(6,*) ’ ENTER: NU2’ 

READ(5,*)NU2 
NU2A-NU2 
TIMEAB*9 . 0E20 

WRITE(6,*) ’ ENTER: IPANEL (O)-RINGS (l)-PANELS’ 

READ (5 , * ) IPANEL 
ISTART-0 

WRITE (6,*) ’ ENTER: LINEUP O-NO SEGMENT LINEUP 1-LINEUP’ 

READ (5,*) LINEUP 

WRITE(6,*) ’ ENTER: MINSECT (MIN # OF THETA SECT. IN SPIRAL)’ 

READ (5 , *)MINSECT 
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WRITE(6,*)’ ENTER: MAXSEGS (MAX # OF RINGS OR PANELS)’ 
READ(5,*)MAXSEGS 
DISTIN-9.0E5 
DISTMX*DISTIN**2 
C 

NSMALL-NRING-1 

LASTSTEP=1 

C 

C — Maximum circulation 
C 

IF (IAXI .EQ. 1) THEN I...3 
GCONST-4.0/PI 
ELSE ! . . .3 
GC0NST*2 . 0 
END IF I . . .3 
C 

C--Compute the initial location oi the tip core 
C 

XO-FRAC 

A0-0.520*X0 

B0-0.826*X0 

C 

C — Set core size 
C 

CCORE- ( 1 . O/NSMALL) /4 . 0 
C 

C — Compute radius, height, and ring circulation for constant spacing 
C 

CALL SETA 

TIME(l) -( (2 . 0*PI**2) / (9 . 0*AE2) ) *X0** (1.5) 

TIMENOW-TIME(l) 

DR* ( 1 . 0-FRAC) /NSMALL 
IF (ISTART .EQ. 0) THEN 
EXPAND* ( 1 . 0- AO) / ( 1 . O-FRAC) 

ELSE 

EXP AND* 1.0 
END IF 

DO 10 1*1 .NSMALL 
HEIGHT (I) *0.0 
CORE(I)«CCORE 
R1-(I-1)*DR 
R2«I*DR 

RADIUS ( I ) -EXP AND+GETRAD ( ILOC , R1 , R2) 

IF (IWING .EQ. 0) THEN I...6 
CIRC1— GCONST* (1 . 0-SqRT(l . 0-Rl**2) ) 

CIRC2— GCONST* (1 . 0-SQRT(l . 0-R2**2) ) 

ELSE ! . . .6 
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CIRC1— GCONST* (1.0- (R1**IWING) *SQRT ( 1 . 0-Rl**2) ) 

CIRC2— GCONST* (1 .0- (R2**IWING) *SQRT(1 .0-R2**2) ) 

END IF !...6 

GAMMA (I) “CIRC2-CIRC1 

10 CONTINUE 
C 

C--Initialize the arc length array, AXLSEG(I) is the total arc length 
C to the end of segment I. (Note: segment I is the segment which 
C contains ring or panel I) 

C 

AXLSEG (0) -0.0 
DO 11 1*1 .NSMALL-l 
AXLSEG ( I ) * AXLSEG ( I - 1 ) 

$ +SQRT( (RADIUS(I+1) -RADIUS(I) )**2 

$ +(HEIGHT(I+1) -HEIGHT(I) )**2 ) 

11 CONTINUE 

AXLSEG (NSMALL) -AXLSEG (NSMALL- 1 ) 

$ +SQRT( (RADIUS (NSMALL) -RADIUS (NSMALL- 1))**2 

$ + (HEIGHT (NSMALL) -HEIGHT (NSMALL-l) ) **2) 

-Initialize special tip vortex (Index-NRING) 

IF (FRAC .EQ. 0.0) THEN I...10 

-This section lor no circulation in tip vortex initially 

HEIGHT (NRING) *0.0 
RADIUS (NRING) *1.0 
GAMMA (NRING) -0.0 
CORE (NRING) -0.0001 

VOLTIP-2 . 0*PI*PI*RADIUS (NRING) *CORE (NRING) **2 
ELSE 1...10 

IF (IAXI .EQ. 1) THEN !...ll tip parameters lor axi-llow 

HEIGHT (NRING) -BO 
REDG-1 . O-FRAC 
GAMMA (NRING) -(- GCONST) - 

$ (-GCONST* (1 . 0- (REDG** IWING) *SqRT(l . 0-REDG**2) ) ) 

RADIUS (NRING) -1 . O-AO 


Now ligure total K.E. and get core size 


SUM 1-0.0 
SUM2-0.0 
DO 82 1-1, NRING 
W-RADIUS(I) 
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Z-HEIGHT(I) 

DO 82 J-l, NRING 

IF (I .NE. J) THEN ! ... 12 
WP-RADIUS(J) 

ZP-HEIGHT( J) 

R1-SQRT( (Z-ZP) **2+ (W-WP) ** 2 ) 

R2=SQRT( (Z-ZP) **2+ (W+WP) **2) 

LAMBDA* (R2-R1 ) / (R2+R1 ) 

CALL ELLIP (LAMBDA .K.E) 

SUM1-SUM1+GAMMA ( I ) *GAMMA ( J) * (R1+R2) * (K-E) /2 . 0 
ENDIF ! ... 12 

82 CONTINUE 
C 

IF (IWING .EQ. 0) THEN ! ... 13 
T0TKE=4. 0/3.0 

ELSE IF (IWING .EQ. 2) THEN '....13 
TOTKE-O.41 

ELSE IF (IWING .EQ. 4) THEN ! . . . 13 
T0TKE-0.245 
ENDIF ! ... 13 
DO 83 I-l.NSMALL 
W*RADIUS (I) 

83 SUM2-SUM2+GAMMA (I)**2*W* ( ALOG (8 . 0*W/C0RE ( I ) ) - 1 . 75 ) /2 . 0 
TNBAR=T0TKE-SUM2-SUM1 

XTEMP-EXP (2 . 0*TNBAR/ (GAMMA(NRING) **2*RADIUS (NRING) ) +1 . 75) 
CORE (NRING) *8 . 0*RADIUS (NRING) /XTEMP 


-Get volume of tip core 

V0LTIP*2 .0*PI*PI*RADIUS (NRING) *C0RE(NRING)**2 

ELSE I...11 tip parameters lor 2-D How 

HEIGHT (NRING) -BO 
REDG-1 .O-FRAC 
GAMMA (NRING) - (-GCONST) - 

$ ( -GCONST* (1.0- (REDG** IWING) *SQRT ( 1 . 0-REDG**2 ) ) ) 

RADIUS (NRING) -1 . O-AO 
CORE (NRING) -0.0 
VOLTIP-O . 0 
ENDIF ! . . .11 

ENDIF ! . . .10 


--Center ring 

IF (IWING .EQ. 0) THEN ! . . . 15 
GMAX— GCONST 
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ELSE ! ... 15 
GMAX-0.0 
END IF ! ... 15 
RAD IUS (0) *0.0 
HEIGHT (0) *0.0 


--I1 panel case, then calculate some more initial parameters 

IF (IPANEL .Eq. 1) THEN ! . . . 16 
DO 1920 1*1 .NSMALL-l 

SL— SQRT( (RADIUS (I) -RADIUS (I- 1) )**2 
$ +(HEIGHT(I) -HEIGHT ( I- 1) )**2 ) 

SR- SqRT( (RADIUS (I) -RADIUS (I+l))**2 
$ +(HEIGHT(I) -HEIGHT(I+1) )**2 ) 

Al- ( (SL*SR) / (SL-SR) ) * ( (RADIUS ( 1+1 ) -RADIUS (I) ) /SR**2 
$ -(RADIUS(I-l) -RADIUS(I) )/SL**2) 

Bl-( (SL*SR)/ (SL-SR) )*((HEIGHT(I+1) -HEIGHT(I) )/SR**2 
$ -(HEIGHT(I-l) -HEIGHT(I) )/SL**2) 

DRPAN(I)*A1 
DZPAN(I)*B1 
1920 CONTINUE 

A2- ( (RAD IUS ( 1+ 1 ) -RADIUS ( I) ) /SR- (RAD IUS ( I - 1 ) -RADIUS ( I ) ) /SL) 

$ /(SR-SL) 

B2*( (HEIGHT(I+1) -HEIGHT(I))/SR- (HEIGHT(I-l) -HEIGHT(I) )/SL) 

$ /(SR-SL) 

DRPAN (NSMALL) *A1+A2*SR 
DZPAN (NSMALL) *B1+B2*SR 
DO 717 1-1, NSMALL 

WPAN(I)*(AXLSEG(I)-AXLSEG(I-l))/2.0 IWPAN is panel half-width 
717 CONTINUE 
END IF ! ... 16 


DO 8315 I*0,NRING 

8315 HEIGHT(I)*HEIGHT(I) -1 .0*TIME(1) 
GMTO*GAMMA(NRING) 


-Write out the input data 


WRITE (IWR,*)' 

NRING 

= ’ , NRING 

WRITE (IWR,*)' 

FRAC 

*’ ,FRAC 

WRITE (IWR,*)* 

INTKOD 

-* , INTKOD 

WRITE (IWR,*)* 

DXMAX 

** , DXMAX 

WRITE(IWR,*) * 

NSTEP 

*’ .NSTEP 

WRITEUWR,*) * 

THMAX 

= ’ , THMAX 

WRITE (IWR,*)' 

IWING 

*' .IWING 

WRITE (IWR,*) * 

IREDIS 

-* , IREDIS 

WRITE (IWR,*)* 

NU2 

-* ,NU2 
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WRITE (IWR,*) ' IPANEL =', IPANEL 
WRITE (IWR,*)’ LINEUP -'.LINEUP 
WRITE (IWR.*) ' MAXSEGS-' .MAXSEGS 
WRITE ( IWR , * ) ' MINSECT* * . MINSECT 
WRITE (IWR,*) * IAXI -'.IAXI 
WRITE (2) INEW 

WRITE (2) NRING , CCORE . IDIST . FRAC . INTKOD . DXMAX . NSTEP . THMAX . 

$ LINEUP . DISTIN . MAXSEGS . MINSECT , IAXI , IPANEL 

CALL VELOCY 
WRITE(IWR,*) * 

WRITE(IWR.*) ' INITIAL VALUES 1 

WRITE(IWR.*) ' 

WRITE (IWR, 3243) 

3243 FORMAT ( * I RADIUS HEIGHT GAMMA 

$ ’ CORE UVEL WVEL AXLSEG 

$ ' DRPAN DZPAN ’ ) 

WRITE (IWR, 3244) (I,RADIUS(I) ,HEIGHT(I) ,GAMMA(I) ,CORE(I) . 

$ UVEL(I) ,WVEL(I) .AXLSEG (I) .DRPAN(I) , 

$ DZPAN(I) ,1-0, NRING) 

3244 FORMAT (IX , 13 , 9E14 . 6) 

--Initialize unit 2 plot data 

WRITE(2) NRING, CORE (NRING) .TIME (1) 

WRITE (2) (RADIUS(I) .HEIGHT(I) .-GAMMA(I) ,AXLSEG(I) .1=0, NRING) 
WRITE (2) TIPVAX 
C 

C 

C MAIN TIME LOOP 

C 

c 

DO 99 ISTEP-LASTSTEP , NSTEP+LASTSTEP- 1 
C 

C — Compute velocities at step ISTEP which 

C corresponds to time (ISTEP-1)* M DT" and also update positions 
C 

CALL MOVEM( ISTEP) 

IF (IREDIS .Eq. 1) THEN 
WVTIP-WVEL (NRING) 

UVTIP-UVEL (NRING) 

CTIP-CORE (NRING) 

NIN-NSMALL 

NOUT-NSMALL 

GMM-GMAX - GAMMA (NRING) 

CALL REDISC (RADIUS .HEIGHT , GAMMA , NIN , NOUT , GMM . SUML . DELTAL , 

$ TIME(ISTEP+1) ) 

WVEL (NRING) -WVTIP 
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UVEL (NRING) *UVTIP 
CORE (NRING) *CTIP 
END IF 

IF (IREDIS .EQ. 1) THEN 
WVTIP=WVEL (NRING) 
UVTIP-UVEL (NRING) 
CTIP-CORE (NRING) 

CALL LUMPS (TIMEUSTEP+1)) 
WVEL (NRING) =WVTIP 
UVEL (NRING) =UVTIP 
CORE (NRING) =CTIP 
END IF 


-Now update the core size, since the sheet has stretched 

IF (IAXI .EQ. 1) THEN 
DO 691 11*1 ,N SMALL 
691 C0RE(II)=WPAN(II)/2.0 
END IF 


--Write data to disk (unit 2) lor plotting later 
NU2-NU2C 

IF (TIME(ISTEP+1) .LT. TIMEBC) NU2=NU2B 
IF (TIME(ISTEP+1) .LT. TIMEAB) NU2=NU2A 

IF ( ( ISTEP/NU2) *NU2- ISTEP . Eq . 0 
$ .OR. ISTEP . EQ . NSTEP+LASTSTEP - 1 ) THEN 

WRITE (2) NRING , CORE (NRING) . TIME ( ISTEP+1 ) 

WRITE(2) (RADIUS (I) .HEIGHT(I) , GAMMA(I) ,AXLSEG(I) . I-O, NRING) 
WRITE (2) TIPVAX 
END IF 

NOWSTEP-ISTEP 
99 CONTINUE 


END OF MAIN TIME LOOP 


STOP 

END 

FUNCTION GETRAD(IL0C,R1,R2) 

C 

C--This function determines the location of the vortex ring within 
C a segment, if IL0C*1 then the ring is centered, if IL0C=2 then the 
C ring is positioned such that the total impulse is maintained 


143 



no o 


IF (ILOC .EQ. 1) THEN 
GETRAD-(Rl+R2)/2.0 
ELSE 

Fl-( (SQRT(1 . 0-Rl**2) ) **3- (SQRTCl . 0-R2**2) ) **3) *8 . 0/3 . 0 
F2»4.0*R1**2*SQRT(1.0-R1**2) 

F3»4 . 0*R2**2*SQRT ( 1 . 0-R2**2) 

F4*4.0*(SQRT(1 .0-Rl**2) -SQRT(1 .0-R2**2) ) 

GETRAD*SQRT ( (F1+F2-F3) /F4) 

END IF 

RETURN 

END 


SUBROUTINE MOVEM(ISTEP) 

COMMON /B1 / RADIUS (0:900) .HEIGHT (0:900) ,WVEL (0:900) ,UVEL (0:900) , 

$ GAMMA (0:900), CORE (0 : 900) , TIME (0 : 5000) , TENERG (0 : 5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING , PI , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /RK / RAD0(0:900) ,HEI0(0:900) ,RK(0:900,4) ,HK(0:900,4) 

COMMON /NOW/ TIMENOW , BIGRREF 
COMMON /FI/ ALPLAST , GMTO 
COMMON /Gl/ TIPVAX 
C 

C- -Computes velocities and new positions using either Euler or 
C Runge-Kutta time stepping 
C 

C- -First save original data 
C 

ALPO* ALPLAST 
DO 20 1-0, NRING 
RADO(I)-RADIUS(I) 

HEIO(I)-HEIGHT(I) 

20 CONTINUE 

-Now check lor time step method 

IF (INTKOD .Eq. 0) THEN ! Euler time step 

TIMENOW-TIME(ISTEP) 

CALL VELOCY 
CALL GETDT(DT) 

CALL GETK(l) 

CALL UPDATE(DT.l) 

TIPVAX-WVEL (NRING) 
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ELSE IRunge-Kutta time step 

C 

CALL PULANG(NRING-1, BEFORE) ! angle at last panel center at start 
TIMENOW-TIME(ISTEP) 

CALL VELOCY 
CALL GETDT(DT) 

CALL GETK(l) 

CALL UPDATE (DT/2. 0,1) 

T IPVAX-WVEL (NRING ) 

C 

CALL PULANG (NRING- l.ANOW) !get new angle and compute the change 
ALPLAST-ALPO+ (ANOW-BEFORE) land add it to the angle of the edge 
TIMENOW-TIME ( ISTEP) +DT/2 . 0 
CALL VELOCY 
CALL GETK(2) 

CALL UPDATE (DT/2. 0.2) 

C 

CALL PULANG (NRING- l.ANOW) 

ALPLAST-ALPO+ (ANOW-BEFORE) 

TIMENOW-TIME ( ISTEP ) +DT/2 . 0 
CALL VELOCY 
CALL GETK(3) 

CALL UPDATE (DT, 3) 

C 

CALL PULANG (NRING- l.ANOW) 

ALPLAST * ALPO+ (ANOW-BEFORE) 

TIMENOW-TIME ( ISTEP ) +DT 
CALL VELOCY 
CALL GETK(4) 

C 

C- -Update positions lor Runge-Kutta 
C 

DO 10 1=0 . NRING 

RADIUS(I)-RAD0(I)+(DT/6.0)*(RK(I , l)+2 .0*RK(I ,2)+ 

$ 2.0*RK(I,3)+RK(I,4)) 

HEIGHT ( I ) -HEIO ( I ) + (DT/6 . 0) * (HK ( I . 1 ) +2 . 0*HK ( I . 2) + 

$ 2 . 0*HK (1,3) +HK (1,4)) 

10 CONTINUE 
C 

END IF 
C 

C--Compute the new core radius lor the tip ring and update time 
C 

CORE (NRING) -SQRT (VOLTIP/ (2 . 0*RADIUS (NRING) ) ) /PI 
TIME ( ISTEP+1) -TIME ( ISTEP) +DT 
C 

RETURN 
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END 

SUBROUTINE GETK(N) 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900) ,CORE(0:900) ,TIME(0: 5000) .TENERG (0:5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING , P I , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /RK/ RAD0(0:900) ,HEIO (0:900) ,RK (0:900, 4) ,HK (0:900, 4) 

C 

C- -Computes the 4 K parameters used in the Runge-Kutta routine. 

C The K's are just the velocity components. 

C 

DO 10 1*0, NRING 
RK(I,N)-UVEL(I) 

HK(I,N)-WVEL(I) 

10 CONTINUE 
C 

RETURN 

END 

SUBROUTINE UPDATE (H.N) 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) , WVEL (0 : 900) ,UVEL(0:900) , 

$ GAMMA (0:900), CORE (0:900), TIME (0:5000), TENERG (0:5000), 

$ TIMPUL (0: 5000) 

COMMON /B2/ NRING, PI, DXMAX. INTKOD, NSMALL, VOLTIP. IAXI 
COMMON /RK/ RAD0(0:900) ,HEI0(0:900) ,RK(0:900,4) ,HK(0:900,4) 

COMMON /PANEL/ IPANEL, DRP AN (0 : 900) ,DZP AN (0:900) ,WP AN (0:900) 

C 

C- -Computes new positions of the rings based on the time step H said 
C the velocity contained in RK and HK. 

C 

DO 10 1-0, NRING 

RAD IUS ( I ) *RADO ( I ) +H*RK ( I . N ) 

HEIGHT ( I ) -HEIO ( I ) +H*HK ( I . N) 

10 CONTINUE 
C 

C--Compute the new core radius lor the tip ring 
C 

CORE (NRING) -SQRT (VOLTIP/ (2 . 0*RADIUS(NRING) ) ) /PI 
C 

C--I1 panel case, estimate inclinations 
C 

IF (IPANEL .EQ. 1) THEN 
DO 20 1*1. NSMALL- 1 

SL— SQRT( (RADIUS(I)-RADIUS(I-1))**2 
$ +(HEIGHT(I)-HEIGHT(I-1))**2 ) 

SR- SQRT( (RADIUS (I) -RADIUS (I+l))**2 
$ + (HEIGHT (I) -HEIGHT(I+1))**2 ) 

Al- ( (SL*SR) / (SL-SR) ) * ( (RADIUS ( 1+1) -RADIUS (I) ) /SR**2 
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$ -(RADIUS(I-1)-RADIUS(I))/SL**2) 

Bl* ( (SL*SR) / (SL-SR) ) * ( (HEIGHT (1+1 ) -HEIGHT (I) ) /SR**2 
$ -(HEIGHT(I-1)-HEIGHT(I))/SL**2) 

DRPAN(I)«A1 
DZPAN(I)*B1 
20 CONTINUE 

A2« ( (RADIUS ( 1+1) -RADIUS (I) ) /SR- (RADIUS (I- 1 ) -RADIUS ( I) ) /SL) 

$ /(SR-SL) 

B2-( (HEIGHT (1+1) -HEIGHT (I) )/SR- (HEIGHT (I- 1) -HEIGHT (I) ) /SL) 

$ /(SR-SL) 

DRPAN (NSMALL) =A1+A2*SR 
DZPAN (NSMALL) *B1+B2*SR 
END IF 

RETURN 

END 

SUBROUTINE GETDT(DT) 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA(0:900) ,CORE(0:900) ,TIME(0:5000) ,TENERG(0: 5000) , 

$ TIMPUL (0:6000) 

COMMON /B2/ NRING, PI, DXMAX.INTKOD. NSMALL, VOLTIP.IAXI 
COMMON /RK/ RADO (0:900) ,HEI0(0:900) ,RK(0:900,4) ,HK(0:900,4) 

-Find max velocity and DT. DT* ( ABS (DXMAX) if DXMAX .LT. 0.0) 

IF (DXMAX .LT. 0.0) THEN 
DT-ABS (DXMAX) 

ELSE 

VMAX-UVEL (0) **2+WVEL (0) **2 
DO 40 I»l. NRING 

VRING-UVEL ( I) **2+WVEL ( I) **2 
IF (VRING .GT. VMAX) VMAX-VRING 
40 CONTINUE 

VMAX-SqRT(VMAX) 

DT-DXMAX/VMAX 
END IF 

RETURN 

END 

SUBROUTINE VELOCY 

-Decide between routines lor velocity from rings or panels, for 
axi or for 2-D 

COMMON /B2/ NRING, PI. DXMAX.INTKOD, NSMALL, VOLTIP.IAXI 
COMMON /PANEL/ IPANEL, DRPAN (0:900) ,DZPAN(0:900) ,WPAN(0:900) 
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IF (IAXI .Eq. 1) THEN 
IF (IPANEL .Eq. 0 ) THEN 

CALL INDVEL ! Rings 

ELSE 

CALL PANVEL ! Panels 

END IF 

ELSE !2-D case 

IF (IPANEL .Eq. 0) THEN 

CALL INDVEL2D ! Filaments 

ELSE 

CALL PANVEL2D ! Panels 

ENDIF 
END IF 
C 

RETURN 

END 

SUBROUTINE INDVEL 
C 

C AXISYMMETRIC FILAMENTS (RINGS) 

C--Tliis routine calculates the induced velocity at each ring position 
C due to all the other rings and due to itself . 

C 

COMMON /Bl/ RADIUS(0:900) ,HEIGHT(0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0 : 900) , C0RE(0 : 900) . TIME(0 : 5000) . TENERG (0 : 5000) , 

$ TIMPUL(0 : 5000) 

COMMON /B2/ NRING ,PI .DXMAX , INTKOD .NSMALL , VOLTIP , IAXI 
REAL*4 LAMBDA. K 

DATA A0.A1.A2.A3.A4.B0.B1.B2.B3.B4/ 

$ 1.38629436112. .09666344259, .03590092383, 

$ .03742563713, .01461196212, .5 

$ .12498593597, .06880248676, .03328355346, 

$ .00441787012/ 

DATA CO , Cl , C2 , C3 , C4 , DO ,D1 , D2 , D3 , D4/ 

$ 1.0 . .44326141463, .06260601220, 

$ .04767383546, .01736506451, .0 

$ .24998368310, .09200180037, .04069697526, 

$ .00526449639/ 

C 

C--Loop over all rings, except center ring, get induced velocity at each 
C 

DO 10 1-1, NRING 
W-RADIUS(I) 

Z-HEIGHT(I) 

WVEL(I)«0.0 

UVEL(I)-0.0 

C 

C — Loop over all rings again (except control ring) and determine the 
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contribution of each to the induced velocity at ring I. 

IF (I .EQ. 1) GO TO 25 
DO 20 J-l.I-i 
WP*RADIUS(J) 

ZP=HEIGHT( J) 

Rl*SqRT( (Z-ZP) **2+ (W-WP) **2) 

R2-SQRT( (Z-ZP) **2+ (W+WP) **2) 

LAMBDA* (R2 -R1 ) / (R2+R1 ) 

FK*1 . 0-LAMBDA**2 

K « AO+FK* (Al+FK* (A2+FK* (A3+FK*A4) ) ) 

$ + (BO+FK* (Bl+FK* (B2+FK* (B3+FK*B4) ) ) ) *ALOG( 1 . O/FK) 

E * CO+FK* (Cl+FK* (C2+FK* (C3+FK*C4) ) ) 

$ + (DO+FK* (Dl+FK* (D2+FK* (D3+FK*D4) ) ) ) *ALOG( 1 . O/FK) 

FI* ( (W-WP) /R1+ (W+WP) /R2)*(K-E) 

F2*2 . 0*E*LAMBDA/ (1 . 0-LAMBDA**2) 

F3* (W* (Rl**2-R2**2) +WP* (Rl**2+R2**2) ) / (R1*R2* (R1+R2) ) 
WV-GAMMA(J)*(F1+F2*F3)/(2.0*PI*W) 

WVEL(I)*WVEL(I)+WV 
F4* (Z-ZP) * (R1+R2) / (R1*R2) 

F5*K-E*(1 .0+LAMBDA**2)/(l .0-LAMBDA**2) 

UV*GAMMA(J) *F4*F5/(2 . 0*PI*W) 

UVEL(I)-UVEL(I) -UV 
20 CONTINUE 

25 IF (I .EQ. NRING) GO TO 40 
DO 30 J-I+l, NRING 
WP-RADIUS(J) 

ZP-HEIGHT(J) 

Rl*SqRT( (Z-ZP) **2+ (W-WP) ** 2 ) 

R2*SqRT ( (Z-ZP) **2+ (W+WP) **2) 

LAMBDA- (R2-R1 ) / (R2+R1 ) 

FK*1 . 0-LAMBDA**2 

K * AO+FK* (Al+FK* (A2+FK*(A3+FK*A4))) 

$ + (BO+FK* (Bl+FK* (B2+FK* (B3+FK*B4) ) ) ) *ALOG( 1 . O/FK) 

E - CO+FK* (Cl+FK* (C2+FK*(C3+FK*C4))) 

$ + (DO+FK* (Dl+FK* (D2+FK* (D3+FK*D4) ) ) ) * ALOG ( 1 . O/FK) 

FI* ( (W-WP) /R1+ (W+WP) /R2) * (K-E) 

F2*2 . 0*E*LAMBDA/ ( 1 . 0-LAMBDA**2) 

F3- (W* (Rl**2-R2**2) +WP* (Rl**2+R2**2) ) / (R1*R2* (R1+R2) ) 
WV-GAMMA(J)*(F1+F2*F3)/(2.0*PI*W) 

WVEL(I)*WVEL(I)+WV 

F4-(Z-ZP)*(R1+R2)/(R1*R2) 

F5*K-E*(1 .0+LAMBDA**2)/(l .0-LAMBDA**2) 
UV-GAMMA(J)*F4*F5/(2.0*PI*W) 

UVEL(I) =UVEL(I) -UV 
30 CONTINUE 
40 CONTINUE 
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C — Calculate final velocity including self-induced component 
C 

WV-GAMMA(I) * (ALOG (8 . 0*W/C0RE (1) ) - . 25) / (4 . 0*PI*W) 
WVEL(I)-WVEL(I)+WV 
C 

C — get velocity of higher order terms in tip region, if not tip ring now 
C 

IF (I .NE. NRING) THEN 

CALL DOUBLET (W , Z , VXPRM , VYPRM , VRPRM , VTHPRM) 

WVEL(I)=WVEL(I) +VYPRM 
UVEL ( I ) -UVEL ( I ) + VXPRM 
END IF 
C 

10 CONTINUE 
C 

C — Get velocity contribution at tip due to inner spiral 
C 

CALL CENSPIR(VXP.VYP) 

UVEL(NRING) -UVEL (NRING) +VXP 
WVEL (NRING) -WVEL (NRING ) +VYP 
C 

C--Get induced velocity at center by extrapolating the velocities of 
C points 1 and 2. Fit a parabola (for the velocity) through points 
C 0,1, and 2 by matching the values at points 1 and 2 and setting 
C the slope-0 at r=0. Then the velocity at r-0 (i.e., the center 
C point) is just the constant part of the parabola. 

C 

UVEL(0)-0.0 

WVEL (0) - (WVEL ( 1 ) *RAD IUS ( 2 ) **2 - WVEL (2 ) *RAD IUS ( 1 ) **2 ) 

$ /(RADIUS (2) **2-RADIUS(l)**2) 

C 

C — done 
C 

RETURN 

END 


SUBROUTINE PANVEL 
C 

C AXISYMMETRIC PANELS 

C--This routine calculates the induced velocity at each panel centroid 
C and at the tip ring due to all the other panels and due to itself . 

C 

COMMON /Bl/ RADIUS(0:900),HEIGHT(0:900),WVEL(0:900),UVEL(0:900), 

$ GAMMA (0:900), CORE (0 : 900) , TIME (0 : 5000) , TENERG (0 : 5000) , 

$ TIMPUL (0:5000) 
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COMMON /B2 / NRING , PI , DXMAX , INTKOO , NSMALL , VOLTIP , IAXI 
COMMON /PANEL/ IPANEL, DRPAN (0: 900) ,DZP AN (0:900) ,WPAN (0:900) 
COMMON /VELO/ R.XX.YY.DXT.DYT.WW.G 
COMMON /Dl/ DISTMX 
REAL*4 LAMBDA. K 


— Loop over all panels, except center ring, get induced velocity at each 
DO 10 1-1, NSMALL 

WVEL(I)*0.0 laxial velocity 

UVEL(I)“0.0 ! radial velocity 


--Loop over all panels again and determine the contribution of each to 
the induced velocity at panel I including the self -induced component. 

DO 20 J-l, NSMALL 
R-RADIUS(J) 

G-GAMMA ( J ) / (2 . 0*WPAN ( J) ) 

WW-WPAN(J) 

DXT-DRPAN(J) 

DYT»DZPAN(J) 

XX-RADIUS ( I ) -RAD IUS ( J) 

YY-HEIGHT(I) -HEIGHT(J) 

DISTSQ»XX*XX+YY*YY 
IF (DISTSq .LT. DISTMX) THEN 
IF (I .Eq. J) THEN 
XX - -.1 * WW 
CALL VELIN ( UINA, VINA) 

CALL VELIO ( UIOA, VIOA) 

CALL VELOUT ( UOUTA, VOUTA) 

XX - .1 * WW 

CALL VELIN ( UINB, VINB) 

CALL VELIO ( UIOB, VIOB) 

CALL VELOUT ( UOUTB, VOUTB) 

UIN - (UINA + UINB) / 2. 

VIN - (VINA + VINB) / 2. 

UIO » (UIOA + UIOB) / 2. 

VIO - (VIOA + VIOB) / 2. 

UOUT - (UOUTA + UOUTB) / 2. 

VOUT * (VOUTA + VOUTB) / 2. 

XX-0.0 

ELSE 

CALL VELIN ( UIN, VIN) 

CALL VELIO ( UIO, VIO) 

CALL VELOUT ( UOUT. VOUT) 

END IF 
ELSE 
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UIN-0.0 

VIN-0.0 

UI0-0.0 

VIO-O.O 

CALL VELOUT (UOUT , VOUT) 

END IF 

UVEL ( I ) -UVEL ( I ) - (UIN+UOUT-UIO) 

WVEL (I) -WVEL (I) - (VIN+VOUT-VIO) 

20 CONTINUE 
C 

C--Now get contribution of tip ring 
C 

J-NRING 

W-RADIUS(I) 

Z-HEIGHT(I) 

WP-RADIUS(J) 

ZP-HEIGHT(J) 

R1-SQRT( (Z-ZP) **2+ (W-WP) **2) 

R2-SqRT( (Z-ZP) **2+ (W+WP) **2) 

LAMBDA- (R2-R1 )/ (R2+R1 ) 

CALL ELLIP (LAMBDA, K,E) 

FI- ( (W-WP) /R1+ (W+WP) /R2)*(K-E) 

F2-2 . 0*E*LAMBDA/ (1 . 0-LAMBDA**2) 

F3=(W* (Rl**2-R2**2) +WP* (Rl**2+R2**2) ) / (R1*R2* (R1+R2) ) 
WV-GAMMA ( J) * (F1+F2*F3) / (2 . 0*PI*W) 

WVEL(I)-WVEL(I)+WV 

F4-(Z-ZP)*(R1+R2)/(R1*R2) 

F6-K-E* ( 1 . 0+LAMBDA**2) / ( 1 . 0-LAMBDA**2) 

UV-GAMMA(J) *F4*FB/ (2 . 0*PI*W) 

UVEL(I)-UVEL(I) -UV 
C 

C--get velocity of higher order terms in tip region 
C 

CALL DOUBLET (W , Z , VXPRM , VYPRM , VRPRM , VTHPRM) 

WVEL ( I ) -WVEL ( I ) +VYPRM 
UVEL ( I ) -UVEL ( I ) + VXPRM 
C 

10 CONTINUE 
C 

C--Get induced velocity at center by extrapolating the velocities of 
C points 1 and 2. Fit a parabola (for the velocity) through points 
C 0, 1, and 2 by matching the values at points 1 and 2 and setting 
C the slope-0 at r-0. Then the velocity at r-0 (i.e., the center 
C point) is just the constant part of the parabola. 

C 

UVEL (0) =0.0 

WVEL (0) - (WVEL ( 1 ) *RADIUS (2) **2-WVEL (2) *RADIUS (1 ) **2) 
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$ 


/ (RADIUS (2) **2-RADIUS(l)**2) 


-Get induced velocity at the tip ring 

WVEL(NRING)=0. 0 
UVEL (NRING) =0 . 0 

DO 40 J-l.NSMALL Hoop over all panels 

R-RADIUS(J) 

G-GAMMA ( J) / (2 . 0*WPAN ( J ) ) 

WW-WPAN(J) 

DXT-DRPAN(J) 

DYT-DZPAN(J) 

XX-RADIUS (NRING) -RADIUS(J) 

YY-HEIGHT (NRING) -HEIGHT (J) 

CALL VELIN ( DIN. VIN) 

CALL VELIO ( UIO, VIO) 

CALL VELOUT( DOUT, VODT) 

UVEL(NRING)«UVEL(NRING) - (UIN+UOUT-UIO) 

WVEL (NRING) -WVEL (NRING) - (VIN+VOUT-VIO) 

40 CONTINUE 

-Get tip ring sell-induced component 

WV-GAMMA (NRING) * (AL0G(8 . 0*RADIUS (NRING) /CORE (NRING) ) -0 . 25) 

$ / (4 . 0*PI*RADIUS (NRING) ) 

WVEL (NRING) -WVEL (NRING) +WV 


— Get velocity contribution at tip due to inner spiral 

CALL CENSPIR ( VXP , VYP ) 

UVEL (NRING) -UVEL (NRING) +VXP 
WVEL (NRING) -WVEL (NRING) +VYP 

--done 

RETURN 

END 

SUBROUTINE INDVEL2D 
C 

C This routine is lor the 2-D case, i.e., filaments 
C It calculates the induced velocity at each lilament position 
C due to all the other lilaments and due to itsell. 

C 

COMMON /Bl/ RADIUS (0:900) .HEIGHT (0:900) .WVEL (0:900), UVEL (0:900), 

$ GAMMA (0:900) .CORE (0:900) .TIME (0:5000) .TENERG (0:5000) . 

$ TIMPUL(0 : 5000) 

COMMON /B2/ NRING , PI , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
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COMMON /PANEL/ IPANEL, DRP AN (0 : 900) ,DZP AN (0:900) ,WP AN (0:900) 

COMMON /VELO/ R.XX.YY.DXT.DYT.WW.G 

COMMON /B9/ IWR 

COMMON /SEGS/ AXLSEG (0:900) 

COMMON /NOW/ TIMENOW , BIGRREF 
COMMON /El/ TAUO , BDUM , CDUM , DDUM , TAU2L00P 
COMMON /FI/ ALPLAST , GMTO 
REAL*4 LAMBDA, K 
C 

TPIINV-1 .0/ (2.0*PI) 

C 

C--Loop over all filaments, except center filament, 

C get induced velocity at each 
C 

DO 10 I-l, NRING ! index of control point 

W-RADIUS(I) 

Z “HEIGHT (I) 

WVEL(I)«0.0 laxial component 

UVEL(I)“0.0 (radial component 

IF (IWR .Eq. 25 .AND. I .GE. 68) THEN 
WRITE (IWR,*) ’ LOOP 10: I,W,Z\I,W,Z 

WRITE(IWR,*) ’ LOOP 20: J, GAMMA ( J) , XX, YY.UIN.VIN.UINRF.VINRF' 
END IF 
C 

IF (I .EQ. 1) GO TO 25 

DO 20 J-1,1-1 (points before control point 

XX»W-RADIUS(J) 

YY-Z-HEIGHT(J) 

Q5Q5-XX*XX+YY * YY 
GCON-GAMMA ( J) *TP I INV/q5q5 
UIN-GCON*YY 
VIN“-GCON*XX 
WVEL(I)“WVEL(I)+VIN 
UVEL(I)“UVEL(I)+UIN 
20 CONTINUE 

25 IF (I .EQ. NRING) GO TO 40 

DO 30 J-I+l, NRING (points after control point 

XX-W-RADIUS(J) 

YY-Z-HEIGHT(J) 

Q6Q5-XX*XX+YY*YY 
GCON-GAMMA ( J) *TPIINV/Q5Q5 
UIN-GCON*YY 
VIN— GCON*XX 
WVEL(I)“WVEL(I)+VIN 
UVEL (I) -UVEL ( I ) +UIN 
30 CONTINUE 

40 DO 45 J*l, NRING (contribution from reflection 
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xx-w+rad ius ( j ) 

YY-Z-HEIGHT(J) 

Q6Q5-XX*XX+YY* YY 
GCON— GAMMA (J)*TPIINV/Q5Q5 
UIN-GCON*YY 
VIN--GCON*XX 
WVEL(I)*WVEL(I)+VIN 
UVEL(I)«UVEL(I)+UIN 
45 CONTINUE 

--Get velocity of higher order terms in tip region, if not tip ring 

IF (I .NE. NRING) THEN 

CALL DOUBLET (W , Z , VXPRM , VYPRM , VEPRM , VTHPRM) 

WVEL ( I ) -WVEL ( I ) + VYPRM 
UVEL ( I ) -UVEL ( I ) + VXPHM 
END IF 

10 CONTINUE 

--Get velocity contribution at tip due to inner spiral 

CALL CENSPIR(VXP.VYP) 

UVEL (NRING) -UVEL (NRING) +VXP 
WVEL (NRING) -WVEL (NRING) +VYP 


Get induced velocity at center by extrapolating the velocities of 
points 1 and 2. Fit a parabola (for the velocity) through points 
0, 1, and 2 by matching the values at points 1 and 2 and setting 
the slope-0 at r-0. Then the velocity at r-0 (i.e., the center 
point) is just the constant part of the parabola. 

UVEL (0) -0.0 

WVEL (0) » ( WVEL ( 1 ) *RAD IUS (2) ** 2 - WVEL (2 ) *RAD IUS ( 1 ) * *2 ) 

$ /(RADIUS (2) **2-RADIUS(l)**2) 


--done 

RETURN 

END 

SUBROUTINE PANVEL2D 

This is for 2-D flow, i.e., 2-D panels 
--This routine calculates the induced velocity at each panel centroid 
and at the tip ring due to all the other panels and due to itself . 

COMMON /Bl/ RADIUS(0 : 900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 


155 


ooooo o o o 


) GAMMA (0:900) ,C0EE(0:900) ,TIME(0:5000) .TENERG (0:5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING , PI , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /PANEL/ IPANEL , DRPAN (0 : 900) . DZPAN (0:900), WPAN (0 : 900) 

COMMON /VELO/ R.XX.YY.DXT.DYT.WW.G 
COMMON /Dl/ DISTMX 
REAL+4 LAMBDA ,K 

— Loop over all panels, except center ring, get induced velocity at each 
DO 10 1-1. NSMALL 

WVEL(I)=0.0 laxial velocity 

UVEL(I)*0.0 ! radial velocity 


—Loop over all panels again (except for the control panel which has 
no sell-induced component) and determine the contribution ol each to 
the induced velocity at panel I. 

DO 20 J«1 .NSMALL 
R-RADIUS(J) 

G-GAMMA ( J) / (2 . 0*WPAN ( J) ) 

WW-WPAN(J) 

DXT-DRPAN(J) 

DYT-DZPAN(J) 

XX*RADIUS( I) -RADIUS (J) 

YY-HEIGHT(I) -HEIGHT( J) 

DISTSQ m XX*XX+YY * YY 
IF (I .NE. J) THEN 

IF (DISTSq .LT. DISTMX) THEN 
CALL VEL2D (UIN , VIN) 

ELSE 

CALL VEL2DF (UIN , VIN) 

END IF 
ELSE 
UIN-0.0 
VIN-0.0 
END IF 

DXT—DXT ! contribution Irom reflection 

G— G 

XX-RADIUS (I) +RADIUS ( J) 

DISTSq-XX*XX+YY*YY 
IF (DISTSq .LT. DISTMX) THEN 
CALL VEL2D(UINRF,VINRF) 

ELSE 

CALL VEL2DF (UINRF , VINRF ) 

END IF 

UVEL ( I ) -UVEL (I) -UIN -UINRF 
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WVEL(I)*WVEL(I)-VIN-VINRF 
20 CONTINUE 

-Now get contribution of tip filament 

J-NRING 
WWP AN-1.0 

G-GAMMA ( J ) / ( 2 . 0* WWPAN ) 

WW-WWPAN 

DXT-DRPAN(J) 

DYT-DZPAN(J) 

XX-RADIUS ( I ) -RADIUS ( J) 

YY-HEIGHT ( I ) -HEIGHT ( J ) 

CALL VEL2DF(UIN,VIN) 

DXT--DXT ! contribution from reflection 

G— G 

XX-RAD IUS(I) +RAD IUS(J) 

CALL VEL2DF(UINRF,VINRF) 

WVEL ( I ) -WVEL ( I ) - VIN- VINRF 
UVEL ( I ) -UVEL ( I ) -UIN-UINRF 


— get velocity of higher order terms in tip region 

CALL DOUBLET ( W , Z , VXPRM , VYPRM , VRPRM , VTHPRM) 

WVEL ( I ) -WVEL ( I ) + VYPRM 
UVEL ( I ) -UVEL ( I ) + VXPRM 

10 CONTINUE 

-Get induced velocity at center by extrapolating the velocities of 
points 1 and 2. Fit a parabola (for the velocity) through points 
0, 1, and 2 by matching the values at points 1 and 2 and setting 
the slope-0 at r-0. Then the velocity at r-0 (i.e., the center 
point) is just the constant part of the parabola. 

UVEL (0) =0.0 

WVEL(0)=(WVEL(1)*RADIUS(2)**2-WVEL(2)*RADIUS(1)**2) 

$ /(RADIUS (2) **2-RADIUS(l)**2) 


--Get induced velocity at the tip filament due to all of the panels 

WVEL (NRING) =0.0 
UVEL (NRING) =0.0 

DO 40 J-l.NSMALL Hoop over all panels 

R-RADIUS(J) 

G-GAMMA ( J ) / ( 2 . 0* WPAN ( J ) ) 

WW-WPAN(J) 
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DXT=DRPAN(J) 

DYT-DZPAN(J) 

XX-RADIUS (NRING) -RADIUS( J) 

YY-HEIGHT (NRING) -HEIGHT (J) 

CALL VEL2D(UIN ,VIN) 

DXT— DXT ! contribution from reflection 

G— G 

XX-RADIUS (NRING) +RADIUS ( J) 

CALL VEL2D (UINRF , VINRF ) 

UVEL (NRING) -UVEL (NRING) -UIN-UINRF 
WVEL (NRING) “WVEL (NRING) -VIN-VINRF 
40 CONTINUE 
C 

C — Get the remaining contribution to the induced velocity at the tip 
C from its reflection 
C 

WWPAN-1.0 

G— GAMMA( J) / (2 . 0*WWPAN) 

WW-WWPAN 
DXT— DRPAN(J) 

DYT- DZPAN(J) 

XX=RADIUS(I)+RADIUS(J) 

YY-HEIGHT(I) -HEIGHT( J) 

CALL VEL2DF (UINRF, VINRF) 

WVEL(NRING) -WVEL (NRING) -VINRF 
UVEL(NRING) -UVEL (NRING) -UINRF 

— Get velocity contribution at tip due to inner spiral 

CALL CENSPIR(VXP.VYP) 

UVEL (NRING) -UVEL (NRING) +VXP 
WVEL(NRING) -WVEL (NRING) +VYP 

— done 
C 

RETURN 

END 

SUBROUTINE IMPULS(TIMP) 

C 

C — Computes total impulse of the vortex filament system for 
C either axi or 2-D. (Note: only computing the "x" or "radial" 

C component of the 2-D impulse. 

C 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0 : 900) , CORE (0 : 900) , TIME (0 : 5000) , TENERG (0 : 5000) . 

) TIMPUL (0:5000) 

COMMON /B2/ NRING , PI .DXMAX , INTKOD , NSMALL , VOLTIP . IAXI 
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IF (IAXI .EQ. 1) THEN 
TIMP-0.0 
DO 10 1*1 ,NRING 

10 TIMP*TIMP+PI*GAMMA(I) *RADIUS(I)**2 
ELSE 

TIMP*0.0 

DO 20 1*1 ,NRING 

20 TIMP“TIMP+GAMMA ( I ) *RADIUS ( I ) 

END IF 

RETURN 

END 

SUBROUTINE TOTENG(ENERG) 


--Computes total energy of the vortex ring system 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900) .CORE (0:900) .TIME (0:5000) .TENERG (0:5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING .PI .DXMAX , INTKOD , NSMALL . VOLTIP , IAXI 
REAL *4 LAMBDA.K 

SUM1-0.0 
SUM2-0.0 
DO 20 1*1, NRING 
W-RADIUS(I) 

Z*HEIGHT(I) 

DO 20 J-l, NRING 
IF (I .NE. J) THEN 
WP-RADIUS(J) 

ZP-HEIGHT(J) 

Rl-SqRT ( (Z-ZP) **2+ (W-WP) **2) 

R2-SQRT ( (Z-ZP) **2+ (W+WP) **2) 

LAMBDA- (R2-R1) / (R2+R1) 

CALL ELLIP (LAMBDA. K.E) 

SUM1*SUM1+GAMMA(I) *GAMMA(J) * (R1+R2) *(K-E) /(2 . 0*PI) 

END IF 

20 CONTINUE 
SUM1-SUM1*PI 

DO 30 1*1, NRING 
W-RADIUS(I) 

30 SUM2-SUM2+GAMMA ( I ) **2*W* (ALOG (8 . 0*W/C0RE ( I ) ) - 1 . 75) / (2 . 0*PI ) 

SUM2=SUM2*PI 

I 

ENERG-SUM2+SUM1 
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RETURN 

END 

SUBROUTINE ELLIP (RK , SUMK , SUME) 

C 

C--Computes complete elliptic integrals, function of RK. 

C Uses a four term series approximation. 

C SUMK - K(RK) - FIRST KIND 

C SUME « E(RK) - SECOND KIND 

C 

DATA AO , A1 , A2 , A3 , A4 , BO , B1 , B2 , B3 , B4/ 

$ 1.38629436112, .09666344259, .03590092383, 

$ .03742563713, .01451196212, .5 

$ .12498593697, .06880248576, .03328355346, 

$ .00441787012/ 

DATA C0.C1,C2,C3,C4,D0,D1.D2.D3,D4/ 

$ 1.0 . .44325141463, .06260601220, 

$ .04757383646, .01736506451, .0 

$ .24998368310, .09200180037, .04069697526, 

$ .00526449639/ 

C 

FK«1.0-RK**2 

SUMK - AO+FK* (Al+FK* (A2+FK* (A3+FK*A4) ) ) 

$ + (BO+FK* (Bl+FK* (B2+FK* (B3+FK*B4) ) ) ) *ALOG ( 1 . O/FK) 

SUME - C0+FK*(C1+FK* (C2+FK* (C3+FK*C4) ) ) 

$ ♦ (DO+FK* (Dl+FK* (D2+FK* (D3+FK*D4) ) ) ) *ALOG ( 1 . O/FK) 

C 

RETURN 

END 

SUBROUTINE SETA 

COMMON /B2/ NRIN6 ,PI .DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /E2/ AO.Al , A2 ,A3 .RHOSTAR, A 
COMMON /El/ TAUO , BDUM , CDUM , DDUM , TAU2L00P 
COMMON /FI/ ALPLAST.GMTO 
COMMON /NOW/ TIMENOW . BIGRREF 
C 

ALPLAST*PI/2 . 0 
TAU0-0.04S 
BIGRREF s 0 . 0966063 
TAU2L00P-1.4542 
IF (IAXI .EQ. 0) THEN 
A«SQRT(2.0) 

ELSE 

A*2 . 0*SQRT (2 . 0) /PI 
END IF 
C 

RETURN 
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SUBROUTINE REDISC (R , Z , DGAM . NIN . NOUT . GMAX , SUML , DELTAL , TIMEY) 

C 

C--This does curve fitting using blended parabolas on R, Z, DGAM 
C 

COMMON /B2/ NRING , P I , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /B9/ IWR 

COMMON /Cl / LINEUP, MAXSEGS.MINSECT 
C 

DIMENSION R(0:900) ,Z (0:900) .DGAM (0:900) ,S (0:900) ,QR (0:900) . 

$ qZ (0:900) , AO (0:900) ,A1 (0:900) ,A2(0:900) ,A3(0:900) . 

$ BO (0:900) ,B1 (0:900) ,B2(0:900) ,B3(0:900) ,C0(0:900) , 

$ Cl (0 : 900) , C2 (0 : 900) , C3 (0 : 900) , GMID (0:900), SMID (0 : 900) , 

$ Q (0:900) ,RS(0:8100) ,ZS(0:8100) ,AL(0:8100) , 

$ THSUB (0:8100) ,ARCG (0:900) .GPL (0:900) 

DIMENSION ROLD (0:900) ,ZOLD(0:900) .DGAMOLD (0:900) 

COMMON/STAR/ ICIRC( 10) .FCIRC(IO) ,SSTAR(10) ,GSTAR(10) .NSTAR 
COMMON /PANEL/ IPANEL, DRPAN (0 : 900) ,DZP AN (0:900) ,WP AN (0:900) 
COMMON /SEGS/ AXLSEG (0:900) 

COMMON /FI/ ALPLAST , GMTO 

LOGICAL LINEAR ! LINEAR * .TRUE, for linear circ. fit 

DATA LINEAR/. FALSE./ ! LINEAR - .FALSE, for cubic circ. fit 

DATA NSUB/9/ 

C 

IL00PER»O 

C 

DO 1 1*0. NRING 
ROLD(I)=R(I) 

ZOLD(I)*Z(I) 

1 DGAMOLD(I)-DGAM(I) 

C 

C--Get straight line distances between points, S(I) is the total 
C distance along the straight lines to node I 
C 

S(0)«0.0 
DO 20 1*1. NIN 

20 S(I)»S(I-l)+SqRT((Z(I)-Z(I-l))**2+(R(I)-R(I-l))**2) 

C 

C — Get value of s for the trailing extended segment, this value 
C of S is out at the presumed end of the sheet. 

C 

S (NIN+1) =S (NIN) +0 . 5* (S (NIN) -S (NIN- 1) ) 

C 

C--Get slopes at each interior node, qR(I)*D(R)/DS at node I 
C 
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DO 30 I-l.NIN-l 

S1»(S(I+1)-S(I))/(S(I-1)-S(I)) 

52- 1.0/S1 

53- 1 .0/ (S(I+1) -S(I-l) ) 

QR(I)-((R(I-1)-R(I))*S1-(R(I+1)-R(I))*S2)*S3 
30 QZ(I)®((Z(I-1) -Z(I) )*S1- (Z(I+1) -Z(I) )*S2)*S3 
C 

C--Coei*s for interior segments 
C 

DO 40 1*1 ,NIN-2 
81-1. 0/(S(I+l) -8(D) 

AO(I)*R(I) 

Al(I)*qR(I) 

A2(I)*S1*(-QR(I+1) -2*QR(I)+3*S1*(R(I+1) -R(I) ) ) 
A3(I)«Sl*Sl*(qR(I+l)+QR(I)-2*Sl*(R(I+l)-R(I))) 
BO(I)-Z(I) 

BKD-qz(i) 

B2(I)»Sl*(-qZ(I+l)-2*qZ(I)+3*Sl*(Z(I+l)-Z(I))) 

40 B3(I)»Sl*Sl*(qZ(I+l)+qZ(I)-2*Sl*(Z(I+l)-Z(I))) 

c 

C--Coef’s ior lirst segment 
C 

S1*S(1)-S(0) 

S2*S(2) -S(O) 

83-1.0/ (8(2) -8(1)) 

A0(0) *R(0) 

A1 (0)-( (R(l) -R(0) ) *S2/S1- (R<2) -R(0))*S1/S2)*S3 
A2(0)*( (R(2) -R(0) ) /S2- (R(l) -R(0) ) /SI) *S3 
A3 (0) -0.0 
B0(0)-Z(0) 

B1 (0)- ( (Z (1) -Z (0) ) *S2/S1-(Z(2) -Z (0) ) *S1/S2) *S3 
B2(0)*( (Z(2) -Z(0) )/S2-(Z(l)-Z(0) )/Sl)*S3 
B3(0)-0 .0 
C 

C--Coef ’s lor last segment 
C 

N-NIN 

S1*S(N) -S(N-l) 

S2«S(N-2)-S(N-l) 

S3-1.0/ (S(N-2)-S(N)) 

A0(N-1)*R(N-1) 

A1 (N-l)*( (R(N) -R(N-l) ) *S2/S1“ (R(N-2) -R(N-l) ) *S1/S2) *S3 
A2(N-1)*( (R(N-2) -R(N-l) )/S2- (R(N) -R(N-l) )/Sl) *S3 
A3(N-l)-0.0 
B0(N-1)*Z(N-1) 

B1 (N-l)*( (Z(N) -Z(N-l) )*S2/S1- (Z(N-2) -Z(N-l) )*S1/S2)*S3 
B2(N-1)*.( (Z(N-2) -Z(N-l) )/S2- (Z(N) -Z(N-l) )/Sl)*S3 
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B3(N-1)*0.0 

--Get the angle of the ring with index NSMALL wrt the tip core 
CALL PULANG (NIN, THNIN) 

— Get the "A" coefficient for the Kadin spiral 
RNU-XNU(THNIN) 

IF (ABS(SIN(THNIN)) .GE. 0.5) THEN 

A* (Z (NRING) -Z (NIN) ) / (THNIN** (-RNU) *SIN(THNIN) ) 

ELSE 

A* (R(NRING) -R(NIN) ) / (THNIN** ( -RNU) *C0S (THNIN) ) 

END IF 
C 

C--Get array of R and Z at the subinterval locations, and compute 
C the total length of the curve 
C 

C- -First do the blended parabola segments 
C 

RS(O)*AO(0) 

ZS(0)*B0(0) 

AL(0)*0.0 
DO 50 1*0, NIN- 1 
DO 50 KSUB-l.NSUB 
J*I*NSUB+KSUB 

T*(S(I+1) -S(I))*FLOAT(KSUB)/FLOAT(NSUB) 
RS(J)*A0(I)+T*(A1(I)+T*(A2(I)+T*A3(I))) 

ZS(J)*B0(I)+T*(B1 (I)+T*(B2(I)+T*B3(I)) ) 

50 AL(J)*AL(J-1)+SQRT((RS(J)-RS(J-1))**2+(ZS(J)-ZS(J-1))**2) 

C 

C — The next few sections find the value of theta at the end of the 
C trailing extended segment, i.e., gets the theta corresponding 
C to arc length S(NIN+1) 

C 

C- -First, before doing the numerical integration, need to get a 
C value to use for delta theta 
C 

DS*S(NIN+1)-S(NIN) 

RNU*XNU (THNIN) 

DTG-DS/ (SQRT (RNU**2+THNIN**2) *A*THNIN** ( -RNU- 1.0)) 
DTG*DTG/2.0 
C 

C--Now have delta theta (i.e. , DTG) which should yield several 
C steps in the numerical integration from S(NIN) to S(NIN+1). 

C 

C--Now do the numerical integration and search for theta(NIN+l) 
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ITERS-0 

SEND*S (MIM+1 ) -S (NIN) 

STH*0 . 0 
TH=THNIM 

64 ITERS* ITERS+1 

IF (ITERS .GT. 100) STOP 'ITERS .GT. 100 LOOP 63* 

RNU-XNU(TH) 

R1*R (NRING) - A*COS (TH) *TH** (-RNU) 

Zl-Z (NRING) -A*SIN(TH) *TH**(-RNU) 

RNU«XNU(TH+DTG) 

R2-R (NRING) -A+COS (TH+DTG) * (TH+DTG) ** ( -RNU) 

Z2=Z (NRING) -A*SIN (TH+DTG) * (TH+DTG) ** (-RNU) 

DS=SQRT ( (R2-R1) **2+ (Z2-Z1) **2) 

IF (STH+DS .GE. SEND) GO TO 63 
STH-STH+DS 
TH=* TH+DTG 
GO TO 64 

63 FRAC*(SEND-STH)/DS 
THNIN1-TH+FRAC+DTG 
C 

C--Now can get the sub-interval locations over the extended spiral 
C segment . 

C 

DTH- (THNIN1 -THNIN) /NSUB 
DO 65 KSUB=1 ,NSUB 
THETA-THNIN+KSDB+DTH 
J-NIN*NSDB+KSDB 
RNU-XNU (THETA) 

RS( J)-R(NRING) -A+THETA** (-RNU) *COS (THETA) 

ZS ( J) -Z (NRING) -A+THETA** ( -RNU) *SIN(THETA) 
AL(J) m AL(J-l)+SQRT((RS(J)-RS(J-l))**2+(ZS(J)-ZS(J-l))**2) 

65 CONTINUE 
C 

C — Total length oi curve is SUML 
C 

SUML-AL( (NIN+1) *NSUB) 

C 

689 CONTINUE 
C 

C — Check for lineup ol spiral sheet segments 
C 

IF (LINEUP .Eq. 1) THEN 

GO TO 2010 !line up the segments 

ELSE 

GO TO 2020 !for equal sized segments 

END IF 
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c 

C--This section for equal sized segments 
C 

C--Get new R and Z points at distance DELTAL apart, these are the 
C centroids 
C 

2020 DELTAL*SUML/NOUT 
XLSEG=DELTAL/2.0 
AXLSEG(0)*0.0 
1*1 

DO 60 J*1 ,NSUB*(NIN+1) 

80 IF (XLSEG .LT. AL(J-l) .OR. XLSEG .GT. AL(J)) GO TO 60 
FRAC*(XLSEG-AL(J-1))/(AL(J)-AL(J-1)) 

R(I)*FRAC*(RS(J) -RS( J-l) )+RS( J-l) 

Z(I)*FRAC*(ZS(J) -ZS(J-l) )+ZS( J-l) 

AXLSEG(I)*XLSEG+DELTAL/2.0 {running total of arc length to 

! segment edge . 

WPAN ( I ) *DELTAL/ 2 . 0 {panel half-width 

I-I+l 

XLSEG*XLSEG+DELTAL 
IF (I .GT. NOUT) GO TO 70 
GO TO 80 
60 CONTINUE 
70 CONTINUE 
GO TO 2030 


This part is for lined up segments in the spiral region 


2010 N0UT1*N0UT 

NTOTSUB=NSUB* (NIN+1) 

CALL ANGLES (RS , ZS , NTOTSUB , THSUB ) 

ALPLAST*THSUB (NTOTSUB) {save angle at edge of sheet for DOUBLET 
IF (THSUB (NTOTSUB) .LE. 9.0*PI/2.0) THEN JALFSPIR is the angle 
ALFSPIR*THSUB (NTOTSUB) -3. 0*PI/2.0 !at which the "spiral" 

ELSE {part of the sheet begins 

ALFSPIR*3 . 0*PI !(i.e. t the equal angular 

ENDIF { segments) 

IF (THSUB (NTOTSUB) .LT. 7.0*PI/2.0) GO TO 2020 {check for at 


{least 1.5 loop 

STOTAL*SUML {total length of sheet 

1*1 {find point on sheet where 

3513 IF (THSUB (I) .GT. ALFSPIR) GO TO 3512 !the spiral section is 
1*1+1 {assumed to start, at 

GO TO 3613 ! theta=ALFSPIR 

3512 SFLAT-AL(I) 

SSPIRAL=STOTAL-SFLAT {length of spiral part of sheet 

R3PIQ2*SQRT( (RS ( I) -ROLD (NRING) )**2+(ZS(I) -ZOLD (NRING) ) **2) 
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RAVERAGE=0 . 0 
DO 3511 J“I .NTOTSUB 
3511 RAVERAGE-RAVERAGE 

$ +SQRT( (RS ( J) -ROLD (NRING) ) **2+ (ZS ( J) -ZOLD (NRING) ) **2) 

RAVERAGE-RAVERAGE/ (NTOTSUB - 1+1 ) 

7126 DENOM*l . 0+ (SSPIRAL*R3PIQ2) / (SFLAT*RAVERAGE) 

NSEG=N0UT1 I# of sheet segments to be output 

NSEGFLAT*INT(NSEG/DEN0M+0.5) ! estimate the # of segments on 
NSEGSP IR=NSEG- NSEGFLAT !the flat and spiral parts 

RNLOOP*(THSUB (NTOTSUB) -ALFSPIR)/(2 .0*PI) !# of loops in spiral 

NTHSECT»INT( (FLOAT (NSEGSPIR)/RNLOOP) +0.5) !# of theta sections 
IF (N0UT1 .LT. MAXSEGS) THEN 
IF (NTHSECT .LT. MINSECT) THEN 
ILOOPER* ILOOPER+ 1 

IF (ILOOPER .GT. MAXSEGS) STOP ’ ILOOPER .GT. MAXSEGS’ 
N0UT1-N0UT1+1 
GO TO 7126 
END IF 
END IF 

DTHETA=2.0*PI/NTHSECT I theta increment of a theta section 

1*1 'segment counter for spiral region 

ISEG-NSEG+l-I ! segment index, we are going backwards 

THEDGE*THSUB (NTOTSUB) -DTHETA ledge of current segment (small theta) 
IF (THEDGE .LT. 7.0*PI/2.0) GO TO 2020 ! check for at least 1.5 loop 

THMID*THEDGE+DTHETA/2 . 0 I theta for centroid of current segment 

AXLSEG ( ISEG) =SUML 
DO 3060 J*NT0TSUB,1 ,-l 

3080 IF (THMID .GT. THSUB(J) .OR. THMID .LT. THSUB(J-l)) GO TO 3050 
FRAC- (THMID -THSUB( J-l) )/(THSUB( J) -THSUB(J-l) ) 

R(ISEG)*FRAC*(RS(J) -RS(J-1))+RS( J-l) 

Z(ISEG)=FRAC*(ZS(J)-ZS( J-1))+ZS( J-l) 

THMID*THMID -DTHETA 

3050 IF (THEDGE .GT. THSUB(J) .OR. THEDGE .LT. THSUB(J-l)) GO TO 3060 
FRAC=(THEDGE-THSUB( J-l) )/(THSUB(J) -THSUB(J-l) ) 

AXLSEG ( I SEG - 1 ) =FRAC* (AL(J)-AL(J-1))+AL(J-1) 

THEDGE*THEDGE-DTHETA 

IF (THEDGE .LT. ALFSPIR) GO TO 3070 Idone with spiral part 

1 * 1+1 

ISEG=NSEG+1-I 
GO TO 3080 
3060 CONTINUE 

3070 CONTINUE II contains number of segments in the spiral region 
NSEGSP IR* I I exact value 

NSEGFLAT*NSEG -NSEGSP IR 

now find segments on flat part of sheet 
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--get segment stretching parameters 
SUMI-0.0 

DO 3761 1*1 .NSEGFLAT-l 
3761 SUMI-SUMI+I 

DSLAST*AXLSEG ( ISEG) - AXLSEG ( ISEG- 1 ) 

HHH- (AXLSEG ( ISEG- 1) -NSEGFLAT+DSLAST) 

$ /(SUMI-NSEGFLAT* (NSEGFLAT-l)) 

DSFIX*DSLAST-HHH* (NSEGFLAT-l) 

C 

AXLSEG (0) *0.0 
1-1 

DS I*DSF IX+HHH* ( I - 1 ) 

XLSEG*DSI/2 . 0 !arc length to filament 

EDGE*DSI !arc length to outboard edge of segment 

DO 3061 J*1 .NTOTSUB 

3081 IF (XLSEG .LT. AL(J-l) .OR. XLSEG .GT. AL(J)) GO TO 3061 
FRAC*(XLSEG-AL(J-1))/(AL(J)-AL(J-1)) 

R(I)*FRAC*(RS(J) -RS( J-1))+RS( J-l) 
Z(I)*FRAC*(ZS(J)-ZS(J-1))+ZS(J-1) 

AXLSEG (1) -EDGE 
I-I+l 

DS I «DSF IX+HHH* ( I - 1 ) 

XLSEG-EDGE+DSI/2.0 

EDGE-EDGE+DSI 

IF (I .GT. NSEGFLAT) GO TO 3071 
GO TO 3081 
3061 CONTINUE 
3071 CONTINUE 

C--set WPAN for all segments 
DO 3110 ISEG-l.NSEG 

WPAN (ISEG) - (AXLSEG ( ISEG) -AXLSEG ( ISEG- 1 ) ) /2 . 0 
3110 CONTINUE 

R (N0UT1+ 1 ) *ROLD (NRING) 

Z (N0UT1+1) -ZOLD (NRING) 

DGAM (N0UT1+ 1 ) -DGAMOLD (NRING) 

N0UT-N0UT1 
NRING-N0UT1+1 
NSMALL-NRING-1 
GO TO 2030 
C 

2030 CONTINUE 
C 

C--H panel case, estimate inclinations (only nec . for R-K) 

C 

IF (IPANEL .EQ. 1) THEN 
DO 920 1*1 .NSMALL-l 
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SL— SQRT((R(I)-R(I-1))**2+(Z(I)-Z(I-1))**2) 

SR- SQRT((R(I)-R(I+1))**2+(Z(I)-Z(I+1))**2) 

A1P- ( (SL*SR) / (SL-SR) ) 

$ *((R(I+1)-R(I))/SR**2-(R(I-1)-R(I))/SL**2) 

B1P- ( (SL*SR) / (SL-SR) ) 

$ *( (Z(I+1) -Z(I))/SR**2- (Z(I-l) -Z(I) )/SL**2) 

DRPAN(I)-A1P 
DZPAN(I)-B1P 
920 CONTINUE 

A2P-((R(I+1) -R(I))/SR- (R(I-l) -R(I) )/SL)/(SR-SL) 
B2P-((Z(I+1) -Z(I))/SR- (Z(I-l) -Z(I) )/SL)/(SR-SL) 

DRPAN (NSMALL) =A1P+A2P*SR 
DZPAN (NSMALL) -B1P+B2P*SR 
END IF 

-Get values of S at the midpoints of the original ring locations 

SMID(0)-0.0 
SMID (NIN) -SUML 
DO 90 I-l.NIN-l 

90 SMID ( I ) - (AL (NSUB* ( 1+ 1 ) ) +AL (NSUB* I ) ) /2 . 0 


-Get values of circulation at the midpoints of the original 
ring locations 

GMID (0) -GMAX 
DO 100 I-l.NIN 

100 GMID(I)-GMID(I-1)+DGAM(I) 

-Curve fit for circulation, test for linear or cubic fit 
IF (LINEAR) THEN 

-The following block is for a linear fit of circulation 

DO 411 I-O.NIN-1 
CO(I)-GMIDU) 

C1(I)»(GMID(I+1)-GMID(I))/(SMID(I+1)-SMID(I)) 

C2(I)=0.0 
411 C3(I)-0.0 

ELSE 

--The following block is for a blended parabola fit of circulation 
--Get slopes at each interior point 
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DO 110 1=1 .NIN-l 

S1=(SMID(I+1)-SMID(I))/(SMID(I-1)-SMID(I)) 

110 q(I)«((GMID(I-l)-GMID(I))*Sl- (GMID(I+1) -GMID(I))/S1) 
$ /(SMID(I+1)-SMID(I-1)) 


— Coel’s lor interior points 

DO 120 I-l.NIN-2 
S1=1.0/(SMID(I+1)-SMID(I)) 

G1=(GMID(I+1) -GMID (I) ) *S1 
COCI)-GMID(I) 

C1(I)*Q(I) 

C2(I)*S1*(-Q(I+1) -2*Q(I)+3*G1) 

120 C3(I)=Sl*Sl*(q(I+l)+q(I)-2*Gl) 

— Coel’s lor lirst segment 

S1*SMID(1) -SMID(O) 

S2=SMID(2) -SMID(O) 

S3=l . 0/ (SMID (2 ) - SMID ( 1 ) ) 

C0(0)=GMID(0) 

Cl (0) »S3* ( (GMID (1) -GMID (0) ) *S2/S1- (GMID (2) -GMID (0) ) *S1/S2) 
C2 (0) -S3* ( (GMID (2) -GMID (0) ) /S2- (GMID ( 1 ) -GMID (0) ) /SI ) 
C3(0)*0.0 


Coel’s lor last segment 

51- SMID (MIN-2) -SMID (NIM-1 ) 

52- SMID (MIN) -SMID (NIN- 1 ) 

S3=l .0/ (SMID(NIN-2) -SMID (MIN) ) 

Gl=GMID(MIM-2) -GMID(NIH-l) 

G2-GMID (MIN) -GMID (NIN-l) 

C0(MIM-1)*GMID(NIN-1) 

C1(HIM-1)=S3*(G2*S1/S2-G1*S2/S1) 

C2 (NIM- 1) =S3* (Gl/Sl -G2/S2) 

C3 (NIN-l) -0.0 

C 

END IF 
C 

651 CONTINUE 
C 

C — Interpolate lor new delta circulation values. 

C XLSEG will now be an array containing the arc length to the 
C edge ol segment I, call it AXLSEG(I) 

C 

1=1 

GPREV-GMAX 
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DO 130 J=1 ,NIN 

150 IF(AXLSEG(I) .LT. SMID(J-l) .OR. AXLSEG(I) .GT. SMID(J))GO TO 130 
T- AXLSEG ( I ) - SMID ( J - 1 ) 

GNEW=C0(J-1)+T*(C1(J-1)+T*(C2(J-1)+T*C3(J-1))) 

DGAM(I)“GNEW-GPREV 

GPREV=GNEW 

I-I+l 

IF (I .GT. NOUT-1) GO TO 140 
GO TO 150 
130 CONTINUE 
140 CONTINUE 

T-SMID (NIN) -SMID (NIN- 1 ) 

GNEW-CO (NIN- 1 ) +T* (Cl (NIN- 1 ) +T* (C2 (NIN- 1 ) +T*C3 (NIN- 1 ) ) ) 

DGAM(NOUT) -GNEW-GPREV 

--Done 

RETURN 

END 


SUBROUTINE ANGLES (RS , ZS , NTOTSUB , THSUB ) 


--Get angles relative to tip lor each subinterval 

COMMON /Bl/ R(0:900) ,Z (0:900) ,WVEL (0:900) ,UVEL (0:900) , 

$ GAMMA (0 : 900) , CORE (0 : 900) . TIME (0 : 5000) . TENERG (0 : 5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING.PI.DXMAX.INTKOD.NSMALL.VOLTIP.IAXI 
DIMENSION RS (0:8100) ,ZS(0:8100) , THSUB (0: 8100) 

REAL*4 MAG 
C 

RT-R(NRING) !tip values of radius 

ZT-Z(NRING) ! and height 

C 

A— 1.0 
B-0.0 

C-RS(O) -RT 
D«ZS(0)-ZT 
DOT»A*C+B*D 
CROSS«A*D-B*C 

MAG-SQRT( (A*A+B*B) * (C*C+D*D) ) 

ARG-DOT/MAG 

IF (ABS(ARG) .GT. 1.0) ARG-SIGN(1 .O.ARG) 

THSUB (0) -SIGN (1 . 0 , CROSS) *ACOS (ARG) 

C 

DO 100 1=1, NTOTSUB 
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A=RS(I-1)-RT 

B=ZS(I-1)-ZT 

ORS(I)-RT 

D*ZS(I)-ZT 

DOT=A*C+B*D 

CROSS=A*D-B*C 

MAG=SQRT( (A*A+B*B) * (C*C+D*D) ) 

ARG-DOT/MAG 

IF (ABS(ARG) .GT. 1.0) ARG*SIGN (1.0, ARG) 

100 THSUB ( I ) “THSUB ( I - 1 ) +SIGN ( 1 . 0 , CROSS ) * ACOS ( ARG) 

RETURN 

END 

SUBROUTINE PULANG (NIN , THNIN) 

COMMON /Bl/ R (0:900) ,Z(0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900) .CORE (0:900) ,TIME(0:5000) .TENERG (0:5000) . 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING . PI . DXMAX . INTKOD . NSMALL . VOLTIP . IAXI 
COMMON /TEMP/ BETA (0:900) 

REALM MAG 

RT-R(NRING) 

ZT-Z (NRING) 

A— 1.0 

B-0.0 

C*R(0) -RT 

D-Z(0)-ZT 

DOT»A*C+B*D 

CROSS=A*D-B*C 

MAG=SQRT ( (A*A+B*B) * (C*C+D*D) ) 

ARG-DOT/MAG 

IF (ABS (ARG) .GT. 1.0) ARG s SIGN (1.0, ARG) 

BETA (0) =SIGN( 1.0, CROSS) *ACOS( ARG) 

DO 100 1=1. NIN 

A«R(I-1)-RT 

B=Z(I-1)-ZT 

C=R(I)-RT 

D=Z(I)-ZT 

DOT=A*C+B*D 

CROSS=A*D-B*C t 

MAG»SQRT( (A*A+B*B) * (C*C+D*D) ) 

ARG-DOT/MAG 

IF (ABS (ARG) .GT. 1.0) ARG=SIGN (1.0, ARG) 

100 BETA(I)=BETA(I-1)+SIGN(1 . 0. CROSS) *ACOS( ARG) 


171 


ooo o o o ooo 


c 

THNIN-BETA(NIN) 

C 

RETURN 

END 

FUNCTION XNU(THETA) 

DATA PI/3.141592654/ 

XNU-2. 0/3.0 

IF (THETA .LT. 2.0*PI) THEN 

XNU*(1 .0/(2 .0*PI-1 .0) ) - (1 .0/ (THETA- 1 .0)) +2. 0/3.0 

END IF 

RETURN 

END 

SUBROUTINE LUMPS (TIMEY) 

C 

C--This subroutine looks at the rolled up part of the sheet and 
C determines if part of it should be truncated and lumped into 
C the tip vortex ring. 

C 

C 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0 : 900) , CORE (0 : 900) , TIME (0 : 5000) , TENERG (0 : 5000) . 

$ TIMPUL (05000) 

COMMON /B2/ NRING ,PI .DXMAX , INTKOD .NSMALL , VOLTIP . IAXI 
COMMON /TH/ THETA (0 : 900) , THMAX 
COMMON /B9/ IWR 
REAL*4 LAMBDA, K 


-- Get segment angles relative to horizintal axis 
CALL GETANG 

--See if any rings exceed the cut-off angle THMAX 


M*0 ! M will contain the # of rings to be lumped 

I-NSMALL-1 ! Highest segment number 
20 IF (THETA(I) .LT. THMAX) GO TO 10 
M=M+1 
I-I-l 
GO TO 20 

Now M contains the # of rings to be lumped into the core 

10 IF (M .Eq. 0) GO TO 999 ! Return 

C 

C--Have M rings to lump, first find new position of tip by conserving 
C impulse . 
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c 

C GTIP * new tip circulation 

C RTIP * new tip radius 

C HTIP * new tip height 

C 

GTIP-0.0 

DO 30 I-NSMALL-M+1 , NRING 
30 GTIP-GTIP+GAMMACI) 

C 

IF CIAXI .EQ. 1) THEN 
RT2-0.0 

DO 40 I-NSMALL-M+1, NRING 

40 RT2-RT2+GAMMA ( I ) *RAD IUS ( I ) * * 2 
RTIP-SQRT (RT2/GTIP) 

HTIP-0.0 

DO 60 I»NSMALL-M+1 . NRING 

60 HTIP«HTIP+GAMMA(I)*HEIGHT(I)*RADIUS(I)**2 
HTIP-HTIP/RT2 

ELSE 

RT2-0.0 

DO 41 I-NSMALL-M+1, NRING 

41 RT2-RT2+GAMMA(I)*RADIUS(I) 

RTIP-RT2/GTIP 

HTIP-0.0 

DO 61 I-NSMALL-M+1 .NRING 

61 HTIP-HTIP+GAMMA (I) ^HEIGHT ( I ) 
HTIP-HTIP/GTIP 

END IF 
CTIP-0.0 
C 

IF (IAXI .Eq.l) THEN 


Only get new core size lor the axi case 
--Get new tip core size by conserving kinetic energy. 
First compute total kinetic energy ol system. 

CALL TOTENG(TKE) 


--Now get the various sums needed to compute TTBAR 


SUM 1-0.0 

DO 00 1-1 .NSMALL-M 

TIBAR-GAMMA ( I ) * *2*RADIUS ( I ) * ( ALOG (8 . 0*RAD IUS(I) /CORE (I)) 
$ -1.75)/2.0 

60 SUM1-SUM1+TIBAR 
C 

SUM2-0.0 
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DO 70 I-l.NSMALL-M 
DO 70 J-l.NSMALL-M 
IF (I .NE. J) THEN 

R1-SQRT( (HEIGHT (I) -HEIGHT(J) )**2+(RADIUS(I) -RADIUS( J) )**2) 
R2-SQRT ( (HEIGHT ( I ) -HEIGHT (J) ) **2+ (RADIUS ( I ) +RADIUS ( J) ) **2) 
LAMBDA- (R2-R1 ) / (R2+R1 ) 

CALL ELLIP (LAMBDA. K,E) 

SUM2*SUM2+GAMMA(I) *GAMMA( J) * (R1+R2) * (K-E) /2 . 0 
END IF 

70 CONTINUE 
SUM3-0.0 

DO 80 I-l.NSMALL-M 

R1-SQRT( (HEIGHT(I) -HTIP) **2+ (RADIUS (I) -RTIP) **2) 
R2-SQRT((HEIGHT(I)-HTIP)**2+ (RADIUS (I) +RTIP)**2) 

LAMBDA- (R2 -R1 ) / (R2+R1 ) 

CALL ELLIP (LAMBDA, K.E) 

SUM3«SUM3+GAMMA(I) *GTIP* (R1+R2) * (K-E) /2 . 0 
80 CONTINUE 

SUM4-SUM3 


-Get energy contribution o 1 tip. 
TTB AR- TKE - SUM 1 - SUM2 - SUM3 - SUM4 


-Now get core size corresponding to this energy. 

CZZZ-EXP (2 . 0*TTBAR/ (GTIP**2*RTIP) +1.75) 
CTIP-8 . O+RTIP/CZZZ 


END IF 


--Put tip variables into main arrays 

RADIUS (NRING) -RTIP 
HEIGHT (NRING) -HTIP 
GAMMA (NRING) -GTIP 
CORE(NRING) -CTIP 
C 

C- -Compute new volume of tip ring 
C 

VOLTIP-2 . 0*PI*PI*RADIUS (NRING) *CORE(NRING ) **2 
C 

C--Now rediscretize remaining portion oi sheet. 

C 

NIN-NSMALL-M 
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NOUT*NSMALL 
GMM-GMAX- GAMMA (NRING) 

CALL REDISC (RADIUS , HEIGHT , GAMMA , NIN . NOUT . GMM , SUML .DELTAL , TIMEY) 

C 

C--done 

C 

999 RETURN 
END 

SUBROUTINE GETANG 
C 

C--This routine computes the angles of the segments relative to 
C the horizontal axis 
C 
C 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900) ,CORE(0:900) ,TIME(0:5000) .TENERG (0:5000) , 

$ TIMPUL (0:5000) 

COMMON /B2/ NRING , PI .DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
C COMMON /B3/ DMIN (0:900) 

COMMON /TH/ THETA (0: 900), THMAX 
DIMENSION DT (0:900) 

C 

C--Get increment for first angle 
C 

DT (0) -ATAN ( (HEIGHT ( 1 ) -HEIGHT (0) ) / (RADIUS ( 1 ) -RADIUS (0) ) ) 

THETA (0) *DT (0) 

C 

C--Get THETA for the rest 
C 

DO 10 I- 1. NSMALL- 1 
XM=RADIUS(I-1) 

XI-RADIUS(I) 

XP-RADIUS(I+1) 

YM-HEIGHT(I-l) 

YI-HEIGHT(I) 

YP*HEIGHT(I+1) 

C 

C--Get sign of cross product of Vi and Vi+1 
C 

SCROSS-SIGN (1 . 0 , (XI-XM) * (YP-YI) - (XP-XI) * (YI-YM) ) 

VMAGS=SQRT( ( (XI-XM) **2+ (YI-YM) **2) * ( (XP-XI) **2+ (YP-YI) **2) ) 

DOT* (XI-XM)* (XP-XI) + (YI-YM)* (YP-YI) 

C 

C--Watch out for numerical inaccuracy in argument of arccos 
C 

DT(I)*0.0 

IF (ABS (DOT/VMAGS) .LE. 1.0) DT(I)=SCROSS*ACOS(DOT/VMAGS) 
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THETA(I)*THETA(I-1)+DT(I) 
10 CONTINUE 


RETURN 

END 

SUBROUTINE DOUBLET ( W , Z , VXPRM , VYPRM , VRPRM . VTHPRM) 

C 

C — Calculates the higher order velocity components in the spiral region 
C 

C W radius of control point 

C Z height of control point 

C VXPRM V sub x prime in the notes, radial velocity component 

C VYPRM V sub y prime in the notes, axial velocity component 

GMTO 


axial velocity component 


strength of tip at start of run 


COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900), CORE (0 : 900) , TIME (0 : 5000) , TENERG (0 : 5000) , 

$ TIMPUL (0:5000) 

COMMON /B2 / NRING , PI .DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 

COMMON /El/ TAUO , BDUM , CDUM , DDUM , TAU2L00P 

COMMON /NOW/ TIMENOW , BIGRREF 

COMMON /E2/ AO.Al , A2 , A3 .RHOSTAR, A 

COMMON /FI/ ALPLAST , GMTO 

COMMON /ZZ/ ZXO.ZYO 

COMMON /B9/ IWR 

DIMENSION B(7, 2) ,D(7, 2) !B(coef., exponent) 


DATA B /. 0053006800 
t .0002020870 

t -1.38567 

t -1.31621 

DATA D /. 0010074100 
| .0000157322 

t -2.37327 

t -2.17250 

BIGRREF “0 . 0966063 
NCOEF-7 


.0013100600 

.0001387070 

-1.36294 

-1.30109 

.0001546680 

.0000101938 

-2.28212 

-2.22138 


.000674976 , 

.000100337 . 

-1.34587 

-1.28252 

.0000656961, 

.0000115224, 

-2.40883 

-2.80713 


.000319442 


-1.33015 


.0000312179, 


-2.36051 


BIGRPRM S BIGRREF* (A*TIME(1) )** (2 . 0/3 . 0) Iconstant in time 

RPRIME s SQRT ( (W-RADIUS (NRING) ) **2+ (Z-HEIGHT (NRING) ) **2) 

R-RPRIME/BIGRPRM 

SINTST* (Z-HEIGHT (NRING) ) /RPRIME 

COSTST= (W-RADIUS (NRING) ) /RPRIME 

THETAST-ATAN2(SINTST . COSTST) 

VRPRM*0.0 
VTHPRM=0 . 0 
DO 10 N-l.NCOEF 
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RRR=GAMMA (NRING) /GMTO ! adjust lor varying tip strength 

TAUA-TAU2L00P+ALPR0T (TIMENOW) * (0 . 60702015) / (2 . 0*PI) 
BN«B(N f l)*TAUA**B(N,2) 

DN»D (N . 1) *TAUA**D (N . 2) 

DDVR- (BN*COS (N* (THETAST+PI/2 . 0-ALPROT (TIMENOW) ) ) 

$ -DN*SIN(N* (THETAST+PI/2 .0-ALPROT (TIMENOW) ) ) ) 

$ * (N/R** (N+l ) ) *RRR* (A** (4 . 0/3 . 0) *TIMENOW** ( 1 . 0/3 . 0) /BIGRPRM) 

DDVTH- (BN*SIN (N* (THETAST+PI/2 . 0-ALPROT(TIMENOW) ) ) 

$ +DN+COS (N* (THETAST+PI/2 . O-ALPROT (TIMENOW) ) ) ) 

$ * (N/R** (N+l ) ) *RRR* (A** (4 . 0/3 . 0) *TIMENOW** ( 1 . 0/3 . 0) /BIGRPRM) 

VRPRM -VRPRM+DDVR 
10 VTHPRM*VTHPRM+DDVTH 

VXPRM«VRPRM*COSTST-VTHPRM*SINTST 

VYPRM-VRPRM*SINTST+VTHPRM*COSTST 

C 

ALF-ALPROT (TIMENOW) 

C 

RETURN 

END 

FUNCTION ALPROT(TIMENOW) 

COMMON /Bl/ RADIUS(0:900) .HEIGHT (0:900) ,WVEL(0:900) ,UVEL(0:900) , 

$ GAMMA (0:900) ,CORE(0:900) ,TIME(0:5000) , TENERG (0 : 5000) , 

$ TIMPUL(0 : 5000) 

COMMON /B2/ NRING , P I , DXMAX , INTKOD , NSMALL , VOLTIP , IAXI 
COMMON /El/ TAUO ,B , C ,D , TAU2L00P 
COMMON /FI/ ALPLAST , GMTO 
C 

ALPROT-ALPLAST-PI/2 . 0 
C 

RETURN 

END 

SUBROUTINE CENSPIR(VXP.VYP) 

C — Calculates the velocity induced at the spiral center due to the 
C inner spiral region 

COMMON /Bl/ RAD IUS (0 : 900) , HEIGHT (0 : 900) , WVEL(0: 900), UVEL (0:900), 

$ GAMMA (0:900) .CORE (0:900) .TIME (0:5000) , TENERG (0 : 5000) . 

$ TIMPUL (05000) 

COMMON /B2/ NRING, PI. DXMAX, INTKOD. NSMALL, VOLTIP. IAXI 
COMMON /El/ TAUO ,B , C ,D , TAU2L00P 
COMMON /E2/ AO , A1 , A2 , A3 , RHOSTAR , A 
COMMON /NOW/ TIMENOW. BIGRREF 
COMMON /ZZ/ ZXO.ZYO 
C 

TAUOLOOP-O. 24016 
BGROLOOP-O. 32094 

BIGRPRM-BGROLOOP* (A*TIME(1) ) ** (2 . 0/3 . 0) 

CC-A** (4 . 0/3 . 0) *TIMEN0W**(1 . 0/3 . 0) /BIGRPRM 
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ALFR=ALPROT (TIMENOW) 

TAUA-TAUOLOOP+ALFR* (0 . 60702015) / (2 . 0*PI) 

C 

ALF-ALPROT (TIMENOW) 

C 

VXTILD—O . 00029091 *CC*TAUA** ( -2 . 30525) 

VYTILD—0 . 00424515*CC*TAUA** ( - 1 . 45439) 

VXP-VXTILD*COS (ALFR) -VYTILD*SIN (ALFR) 

VYP=VXT ILD * S IN (ALFR) +VYTILD*COS (ALFR) 

C 

RETURN 

END 

SUBROUTINE VELOUT (UO, VO) 

COMMON /VELO/ R, X, Y, DX, DY. W, G 

REAL Ml. LAM, KO, K1 , K2, KLO, KL1 , KL2, K. KP 

UO « 0.0 

VO - 0.0 

PI - 3.141592654 
R1 - SQRT(X**2+Y**2) 

R2 - SqRT((X+2.*R)**2+Y**2) 

LAM - (R2-R1)/(R1+R2) 

Ml « 1. - LAM**2 

IF (Rl.Eq.O) GO TO 730 

IF (LAM.EQ.O) GO TO 730 

DXR1 - X/Rl 

DYR1 - Y/Rl 

DXR2 - (X+2.*R)/R2 

DYR2 - Y/R2 

DXL - (DXR2-DXR1)/(R2+R1) - LAM/(R2+R1)*(DXR1+DXR2) 

DYL - (DYR2-DYR1 ) / (R2+R1 ) - LAM/(R2+R1)*(DYR1+DYR2) 

KO - 1.3862944 

K1 - .1119723 * Ml 

El - .4630151 * Ml 

K2 - .0725296 * Ml* *2 

E2 - .1077812 * Ml**2 

KLO - .5 

KL1 - .1213478 * Ml 

ELI » .2452727 * Ml 

KL2 - .0288729 * Ml**2 

EL2 * .0412496 * Ml**2 

K - KO + K1 + K2 - (KLO + KL1 + KL2)*L0G(M1) 

E - 1. + El + E2 - (ELI + EL2)*L0G(M1) 

KP - E / LAM / Ml - K / LAM 

EP *(E - K) / LAM 

DXR - DXR1 + DXR2 

DYR - DYR1 + DYR2 

RR - R1 + R2 
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730 


C 

C 


740 



C 

750 


UO - G*W/PI/(R+X) * (DYR*(K-E)+RR*(KP-EP)*DYL) 
VO — G*W/PI/(R+X) * (DXR*(K-E)+RR*(KP-EP)*DXL) 
RETURN 
END 

SUBROUTINE VELIN (UIN, VIN) 

COMMON /VELO/ R, X. Y, DX, DY. W, G 
PI - 3.141592654 

COST * DX / SQRT(DX**2+DY**2) 

SINT - DY / SQRT (DX**2+DY**2) 

Q5 - (X**2 + Y**2)/W**2 
Q4 - (Y*COST - X*SINT)/W 
q3 - (X*COST + Y*SINT)/W 
q2 » q5 - 2.*q3 + l. 
qi ■ qs + 2.*q3 + l. 

IF (q4.Eq.O) THEN 
ATN - 0.0 
GO TO 740 
END IF 

ATN - ATAN((l.-q3)/q4)+ATAN((l.+q3)/q4) 
CONTINUE 

U1 — SINT*L0G(q2/qi)+2.*C0ST*ATN 
VI *-C0ST*L0G(q2/qi)-2. *SINT*ATN 
U2 « (X/2./R) * U1 
V2 - (X/2./R) * VI 

+ . 5*W/R *( L0G(q2*qi) 

-q3*L0G(q2/qi) 

+2.*q4*ATN 

-4.) 

UIN — G / 4. / PI * (U1 + U2) 

VIN - G / 4. /PI * (VI + V2) 

RETURN 

END 

SUBROUTINE VELIO (UIO, VIO) 

COMMON /VELO/ R, X, Y, DX, DY, W, G 
PI * 2* ASINCl.) 

COST - DX / SqRT(DX**2+DY**2) 

SINT - DY / SqRT(DX**2+DY**2) 
q5 » (R**2 + Y**2)/W**2 
q4 - (Y*COST + R*SINT)/W 
q3 - (-R*COST + Y*SINT)/W 
q2 - qs - 2*q3 + l. 
qi * qs + 2*q3 + l. 
uio » o.o 
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VIO - 0.0 

c 

Rl « SQRT(X**2+Y**2) 

IF (Rl.Eq.O.O) GO TO 720 

UIO »-G*W/PI * ( 1. + X/R/2.) * Y/R1/R1 

VIO = G*W/PI * (( 1. + X/R/2.) * X/R1/R1 + . 5/R*L0G(Rl/8/R) ) 

720 VIO - G*W/PI * LOGC 8*R/W ) /2 /R 
* + VIO 

RETURN 
C 

END 

SUBROUTINE VEL2D (UIN, VIN) 

C 

COMMON /VELO/ R, X, Y. DX, DY. W. G 
PI - 2*ASIN(1 . ) 

COST - DX / SQRT(DX**2+DY**2) 

SINT « DY / SqRT (DX**2+DY**2) 

Q5 ■ (X**2 + Y**2)/W**2 
Q4 - (Y*COST - X*SINT)/W 
q3 » (X*COST + Y*SINT)/W 

q2 » qs - 2*q3 + 1. 
qi ■ qs + 2*q3 + 1. 

IF (q4.Eq.O) THEN 
ATN * 0.0 

GO TO 740 

END IF 

ATN - ATAN ( ( 1 . -q3) /q4) +ATAN ( ( 1 . +q3 ) /q4) 

740 CONTINUE 

IF ( (qi . Eq . 0) . OR . (Q2 . Eq . 0) ) THEN 
UIN - 0 

VIN - 0 

GO TO 750 

END IF 

U1 *-SINT*LOG (q2/qi) +2 . *CQST*ATN 
VI «-C0ST*L0G(q2/qi)-2.*SINT*ATN 
UIN - -G / 4./ PI * U1 

VIN - G / 4./ PI * VI 

750 RETURN 

END 

SUBROUTINE VEL2DF (UIN. VIN) 

C 

COMMON /VELO/ R. X. Y, DX. DY. W. GG 

G»GG*2.0*W 

PI » 3.141592654 

qs « SqRT(X**2 + Y**2) 

IF (qS.Eq.O) GO TO 750 
COST - x / qs 
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SINT * Y / Q5 
U1 « SINT/q5 
VI » C0ST/q5 
UIN —G/2./PI * U1 
VIN * G/2./PI * VI 
750 RETURN 
END 
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